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I.  General  Introduction 

The  purpose  of  this  report  is  to  describe  the  results  of  a  survey  on 
the  techniques  that  are  available  for  studying  the  effect  of  random  parameters 
on  the  response  characteristics  for  linear  dynamical  systems.  We  are,  there¬ 
fore,  concerned  with  the  properties  of  the  solutions  of  ordinary  linear 
differential  equations  with  random  parameters. 

There  are  two  classes  of  random  parameters  that  we  distinguish  between 
in  this  report.  They  are  time  independent  (constant)  random  parameters, 
and  time  varying  random  parameters.  In  both  cases,  structural  applications 
are  of  interest. 

Before,  we  outline  the  content  of  this  report,  we  briefly  consider  the 
history  of  the  general  problem  of  parameter  uncertainties  and  parameter 
fluctuations  to  see  why  they  cannot  be  overlooked. 

Very  early  in  engineering  design,  factors  of  safety  were  introduced  to 
account  for  our  lack  of  precise  knowledge  of  the  structure  and  its  loads, 
factors  of  safety  quantified  the  fact  that  loads,  material  behavior, 
structural  element  properties,  etc.,  could  not  be  accurately  estimated. 

Indeed  Roebling  and  his  engineers  in  1880's,  by  means  of  very  careful 
computations,  estimated  the  factor  of  safety  of  the  then  new  Brooklyn 
Bridge  to  be  above  rive.  Later  estimates  in  1944  placed  this  figure  at 
four.  These  factors  essentially  tell  us  that  column  misalignment, 
residual  stress  due  to  manufacturing  errors,  reduction  of  working  area  due 
to  corrosion,  member  weight,  joint  behavior,  etc.,  etc.,  would  have  to  be 
large  indeed  before  the  integrity  of  the  structure  could  be  put  into 
question.  In  other  words,  parameter  variability  would  have  to  be  beyond 
all  reasonable  bounds  before  the  bridge  could  fail.  Roebling  and  his 
engineers  were  correct  in  their  estimation  of  the  effect  of  parameter 
variability  on  the  safety  of  their  bridge.  However,  in  the  first  Quebec 
Railway  Bridge  (1907)  errors  in  estimation  of  member  weight,  in  the  effect 
of  column  misalignment  and  in  the  behavior  of  new  types  of  columns,  lead 
to  failure  because  the  effect  of  parameter  variability  exceeded  the  bounds 
that  the  factors  of  safety  could  absorb. 


Along  the  same  lines,  consider  the  buckling  under  axial  loading  of 
thin  cylindrical  and  conical  shells.  The  buckling  load  of  such  shells 
can  be  calculated.  However,  in  1950' s,  1960's  extensive  tests  of  such 
shells  revealed  that  the  actual  buckling  loads  were  significantly  lower 
than  calculated.  This  discrepancy  between  calculated  and  actual  was 
traced  to  the  fact  that  there  were  random  deviations  from  the  regular 
geometrical  shape  assumed.  These  are  reflected  in  the  fact  that  the  PDE's 
whose  solution  produces  the  buckling  load  contains  random  variations,  in 
its  parameters. 


Randomly  time  varying  parameters  occur  to  a  great  extent  as  a  result 
of  environmental  fluctuations  that  effect  the  system.  The  vibrations  in 
aerospace  vehicles  due  to  atmospheric  turbulence  is  one  prime  example. 

This  is  reflected  in  random  load  variations  on  the  structure  which  are  of 
particular  significance,  for  example,  on  rotating  lifting  surfaces  such  as 
helicopter  blades.  There  has  been  a  great  deal  of  attention  to  this 
problem  due  to  the  critical  nature  of  the  stability  and  safety  of 
helicopters. 


Liquid  sloshing  in  the  tanks  which  are  undergoing  vertical  excitations, 
also  is  a  problem  that  was  actively  studied  for  the  stability  of  the  initial 
atmospheric  stages  of  the  lifting  of  large  rockets.  In  both  of  these 
problems,  randomly  fluctuating  parameters  are  present  in  the  analytical 
equations  that  model  the  response  characteristics  of  these  components.  In 
general,  inverted  beams,  pendulums  as  well  as  aerodynamic  panels  subjected 
to  random  end  loads,  will  be  described  by  models  with  randomly  varying 
parameters. 


Moreover,  for  the  control  of  such  systems  with  uncertain  parameters, 
it  is  necessary  to  be  able  to  characterize  the  response  of  the  controlled 
system,  in  order  to  determine  the  accuracy  required  to  achieve  prescribed 
control  accuracies. 


Surveying  the  techniques  and  results  in  this  overall  class  of  problems, 
it  appears  that  there  exists  a  natural  distinction  between  the  case  of 
random  time  dependent  (fluctuating)  parameters  and  the  case  of  the  random 
time  independent  (constant)  parameters.  We  shall,  therefore,  distinguish 
between  these  two  cases. 


We  shall  carry  on  this  distinction  in  the  next  sections  of  the 
Introduction,  as  well  as  the  general  survey  of  results. 


s 


I. 1  Introduction  -  Time  Varying  Structural  Equations 

Time  varying  models  of  engineering  systems  occur  for  a  number  of  geo¬ 
metrical,  environmental  as  well  as  analytical  reasons.  They  occur  for 
geometrical  reasons  as  a  result  of  the  location  and  directions  that  external 
excitations  impact  upon  a  system,  while  the  system's  physical  parameters  (damping, 
stiffness,  etc.)  are  assumed  to  remain  fixed.  They  occur  as  a  result  of  environ¬ 
mental  properties  due  to  chemical  effects,  thermal  effects,  and  radiation  effects 
that  are  reflected  in  varying  physical  parameters  for  the  components  of  the 
system.  These  may  be  of  a  periodic  random  nature  or  of  a  monotonic  random  nature 
due  to  ageing  of  the  components  in  general.  Finally,  they  occur  in  studies  of 
non-linear  systems.  In  particular,  if  one  wishes  to  study  the  linearized 
equations  of  small  oscillations  about  some  specific  (non-equilibrium)  system 
response,  time  varying  coefficients  will  be  present. 

In  the  geometrical  case,  those  systems  which  are  subjected  to  base 
excitations,  such  as  pendulums  and  missiles;  to  end  loadings,  such  as  beams 
with  various  supports,  and  finally  to  boundary  edge  loadings,  such  as  plates  and 
shells,  will  be  described  by  differential  equations  with  time  varying  coefficients. 

These  cases  are  illustrated  in  the  following  figure 


J_V”< 


P  +n(t) 
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Figure  1 
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These  cases  lead  to  differential  equations  of  the  following  form, 

[1.1.  1.2,  1.3]* 


Pendulum 

(a)  0  +  c9  +  (k+n(t)9=0,  (k  <  0  for  inverted  pendulum) 

Single  Supported  Beam 

(b)  w  +  2gwr  +  (p  +n(t))w  +  w  =0 

tt  c  o  xx  xx  xx 

w(o,t)=w  (o,t)=0,  w(l, t)  =  w  (l,t)  =  0 

XX  xx 

Infinite  panel  of  unit  width  in  supersonic  flow 

(c)  w  +  2gw,.  +  Mw  +  (p+h(t))w  +  w  =0 

Ct  t  x  xto  xx  xxxx 

(boundary  conditions  same  as  above) 


(1.1) 


Using  modal  expansions,  (b)  and  (c)  lead  to  equations  of  the  form  (a). 

For  the  general  linear  structure  that  we  shall  be  concerned  with  in  this 
development,  the  model  will  be  assumed  to  be  of  the  form 


M  x+C(t)  x+K(t)x  =  ?(t) 


(1.2) 


The  response  vector  x  ,  will  be  n-dimensional  for  the  n-mass  structure, 
the  nxn  mass  matrix  M,  will  always  assumed  to  be  known  and  most  often  fixed 
(although  this  is  not  necessary).  The  nxn  damping  matrix  C(t)  and  stiffness 
matrix  K(t)  will,  in  general,  contain  randomly  varying  elements.  The  n-vector 
f  represents  external  excitations  that  may  be  random.  The  external  excitation 
vector  f  does  not  pose  any  analytical  difficulties  and  may  be  treated  as  the 
non-homogeneous  part  of  any  linear  differential  equation.  It  is  the  randomly 
fluctuating  coefficients  in  the  matrices  C(t),  K(t)  that  generate  the  dif¬ 
ficulties.  Time  varying  systems  are  difficult  to  analyze  quantitatively  even 
in  the  deterministic  setting.  Thus,  the  random  setting  will  be  at  least  as 
difficult.  Naturally,  all  will  depend  upon  the  assumed  properties  of  these 
random  coefficients. 

*  References  in  Sections  I-V  are  given  after  Section  V. 
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If  these  random  coefficients  are  due  to  fluctuating  external  loadings, 
they  will  have  frequency  spectra  (power  spectral  densities  in  the  second  order 
stationary  case)  that  will  cover  a  band  containing  low  frequency  as  well  as 
high  frequency  components.  If  the  bandwidth  is  large  and  relatively  constant, 
then  in  many  cases  it  makes  sense  to  model  the  coefficient  fluctuations  as 
Gaussian  white  noise.  If,  on  the  other  hand,  there  are  definite  peaks  in  the 
frequency  spectrum  and  the  larger  frequency  components  are  less  pronounced 
(essentially  band-limited)  then  Gaussian  white  noise  is  not  a  suitable  model. 

In  this  case,  so-called,  physical  noise  is  the  proper  model  for  the  coefficient 
fluctuations.  Finally,  when  the  coefficient  fluctuations  are  due  to  environ¬ 
mental  effects  (thermal,  chemical,  radiation,  or  ageing)  it  is  natural  to  model 
the  fluctuations  as  slowly  varying  (i.e.,  narrow  band,  low  frequency).  Small 
parameters  can  be  applied  via  approximation  schemes.  In  a  related  situation 

if  there  are  small  random  fluctuations  about  a  nominal  value,  again  small  parameters 
can  be  applied  via  approximation  schemes.  However,  in  this  last  case,  we  must 
be  careful.  For  example  for  the  simple  undamped  oscillator  [1.4] 


+  en(t) )x  =  0, 


(1.3) 


if  n(t)  is  the  gaussian  white  noise  then  no  matter  how  small  e>0,  the 
second  moments,  as  well  as  the  sample  solutions,  will  become  unbounded.  Thus, 
even  though  random  fluctuations  are  small  in  their  variances  the  system  can 
still  become  unstable.  This  important  point  cannot  be  overstressed  in  analyzing 
systems  with  random  coefficients.  Therefore  careful  attention  must  be  paid 
always  to  meaningful  approximations,  especially  when  small  parameters  are 
present . 

The  importance  of  the  white  noise  assumption  for  the  coefficient 
fluctuations  is  that  the  solution  process  is  Markov  allowing  us  to  use  the 
many  tools  available.  In  particular,  the  Fokker-Planck  equation,  the  generator 
of  the  associated  diffusion  process  and  finally,  the  related  Ito  differential 
formula  can  be  applied. 
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The  statistical  moments  such  as  means,  variances  and  covariances  can  be 
obtained  explicitly  for  these  linear  models.  In  general,  however,  probability 
densities  cannot  be  obtained  explicitly. 

For  physical  noise  coefficient  processes,  the  story  is  quite  different. 
Here,  we  cannot  even  obtain  the  moments  explicitly,  unless  very  specific 
simplifying  assumptions  are  made.  For  example,  if  the  system  coefficient 
matrices  generate  a  Lie  algebra  [1.5]  then  the  solution  moments  can  be  obtained 
explicitly.  If  the  physical  noise  is  almost  a  white  noise,  an  associated  system 
can  be  studied  with  white  noise  coefficients,  having  statistical  properties  that 
are  similar  to  the  original  system  [1.6],  [1.7],  [1.8].  Under  assumptions 
of  small  parameters,  an  associated  Markov  process  can  be  obtained,  which  also 
will  yield  valid  approximations  to  the  statistical  properties  of  the  original 
system.  [1.9],  [1.10],  [1.11]. 

Finally,  in  lieu  of  all  of  these, successive  approximations  must  be 
applied  [1.12-1.18]  simply  based  upon  the  assumed  physical  noise  statistics. 

The  essential  reason  that  the  physical  noise  coefficient  case  cannot  yield 
the  exact  statistical  properties  of  the  response,  is  simply  that  the  co¬ 
efficient  process  at  any  given  time  is  correlated  to  the  response  at  that  time. 
To  illustrate  this  fact,  let  us  consider  the  simple  first  order  scalar  equation 


x  +  (a  +  n(t) )x  =  0, 


(1.4) 


Upon  taking  expectations,  we  find 


—  E(x(t)}  +  aE  (xvw)  +  E(n(t)x(t)} 


(1.5) 


It  is  the  term  E{n(t)  x(t) ,  that  is  of  concern  to  us.  Since  n(t)  is 
a  physical  noise,  with  dependence  upon  the  past,  we  cannot  simplify  this  term 
any  further.  On  the  other  hand,  if  n(t)  is  the  gaussian  white  noise,  then 
n(t ) ,  x ( t )  are  independent  random  variables  for  any  t,  which  allows  us  to  write 


E{n(t)x(t)}  =  Ein(t)]Ex(t) ■ 


(1.6) 


^  j  ^ ■  *  ,»  *  ,  ' 


In  view  of  (1.6),  we  can  Chen  write  (1.5)  as 


a  m(t)  +  (a +  E(n(t)} )m(t)  =  0, 


(1.7) 


which  is  a  simple  linear  scalar  equation  for  the  mean  response, 
m(t)  =  E  {x(t)  }. 

Naturally,  this  extends  to  higher  order  structural  system  equations  as 
well.  It  is  interesting  to  note  at  this  time,  that  the  many  recursive  schemes 
for  obtaining  the  approximations  to  the  moments,  or  probability  densities,  for 
the  physical  noise  case  base  the  initial  approximation  on  the  assumption  of 
independence  of  the  coefficient  process  with  the  response  process.  At  this  time 
there  does  not  appear  to  be  any  general  technique  available  to  obtain  even  the 
exact  moments  for  the  linear  system  with  physical  noise  coefficients. 

There  are  other  interesting  properties  that  must  be  taken  into  account  when 
studying  the  statistical  properties  of  the  dynamics  of  systems  with  uncertain  or 
randomly  fluctuating  coefficients.  It  is  possible  that  the  average,  statistical, 
properties  of  the  response  process  may  be  quite  distinct  in  character  from  the 
actual  sample  solutions  themselves.  In  particular,  the  asymptotic  behavior  of 
the  mean  motion  may  be  divergent  or  even  be  undefined  (so-called  finite 
explosion  time)  yet  the  actual,  sample,  behavior  of  the  response  will  remain 
quite  regular.  This  anomalie  occurs  for  systems  with  random  coefficients  since 
we  are  dealing  not  only  with  the  solutions  of  differential  equations,  but  with 
their  ensemble  averages  as  well. 

We  illustrate  this  non-intuitive  behavior  with  two  simple  first  order 
scalar  equations. 

Example  1.1  -  Finite  explostion  time  for  moments. 

Let  ,b,  be  a  zero  mean  gaussian  variable  (constant  in  time),  with  unit 

b2 

1  2 

variance.  Thus,  p(b)  =  -  e  ,  the  probability  density  for  b. 

/27 
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We  study  the  mean  behavior  of  the  solution  to  the  scalar  equation 


x(t) =  b  x(t) ,  x(0) =  1 


(1.8) 


We  can  easily  find  the  sample  solution  to  (1.8)  to  be 

,2 

x  ( t )  =  e  C ,  t  10 


(1.9) 


Clearly,  the  samples  are  exponentially  increasing  as  tf°°,  with  a  rate 

2 

depending  upon  the  specific  value  b  e  [0,°°). 

Now,  what  can  be  said  about  the  behavior  of  the  moments,  m^(t)  =  E(x  (t)}? 


We  merely  have  to  evaluate  the  expectation 


E{xk(t) }  =  E{ekb  C}  =  — 


2  kb  t,, 
e  e  db 


/<kt-  i>db. 


(1.10) 


1  t  h 

We  see  from  (1.10),  quite  obviously,  if  t>— ,  then  the  k  moment 
does  not  exist!  That  is,  it  becomes  infinite.  Thus,  for  t >  all  moments 

are  infinite,  even  though  the  sample  solution  behavior,  from  (1.9)  is  quite 
regular.  Hence,  the  mean  motion  will  not  exist  after  finite  time,  even  though 
the  system  response  is  well  behaved  as  exponentially  increasing  curves. 

We  now  look  at  a  somewhat  more  complicated  case,  with  randomly  varying, 
physical  noise,  coefficients. 


*» 


-10- 


Example  1.2  -  Increasing  Moments  -  Decreasing  Samples 

For  this  example,  we  assume  the  fluctuating  coefficient,  scalar  first 
order  equation, 

x(t) +  (a  +  n(t)  x(t)  =  0, 

where  a  >0  is  a  known  constant  and  the  stochastic  process  {n(t),  t>0}, 
is  gaussian  with  zero  mean  and  covariance  function 


y  (t,s)  =  E{n(t)n(s)}  =  a2e  ^ 
n 


(1.11) 


(1.12) 


2a 

1  +  oj2 

5 

V, 

one. 

The  n-process  is  a  stationary  gaussian  process,  whose  spectral  density 

9a2 

- »  ,  is  absolutely  continuous.  Therefore,  it  is  known  to  be  [1.19]  ergodic. 

1  +  L0 

Thus,  we  can  equate  time  averages  and  ensemble  averages,  with  probability 

one.  t 

]_  f 

In  particular,  lim  —  n(s)ds =  E{n(* ) } =  0  with  probability  one. 
tt<*>  1  o 

The  n-process  is,  in  fact,  the  so-called  Ornstein-Uhlenbeck  process  [1.20] 
whose  sample  functions  are  known  to  be  continuous  functions  with  probability 
one.  (This  is  in  contrast  to  the  random  telegraphic  signal,  a  non-gaussian 
process  whose  covariance  is  given  by  (1.12),  but  whose  samples  are  piecewise 
constant  -  (see  [1.21]  for  a  discussion))  . 

Thus,  we  can  integrate  (1.11)  directly,  to  yield  the  sample  solutions  for 
x(0)  =  1, 


x(t)  =  e 


-at  -  n (s)ds 


(1.13) 


,  i 


x 

a 


|  i 


In  order  to  study  the  asymptotic  behavior  of  the  response,  x(t),  we  form 
the  limit. 


-at-  n(s)ds  \ 

lim  x(t)  =  lim  e  =  lim  e 

£  'fee  £  -f  oo  t^00 


(  r  ) 

a-t j  n(s)dsy  t 


(1.14) 


-11- 


However,  by  the  ergodic  properties  of  the  n-process,  we  have 


lim  (-a  -  — 
t  +« 


(t 

n(s)ds) 

o 


i  ^  1 

-a  -  lim  — 

-foo 


t 

n(s)ds 

o 


(1.15) 


=  -a,  with  probability  one. 


Finally,  since  e  10  as  t+“,  it  follows  that  lim  x(t)  =  0  with 

ft'® 

probability  one.  This  property  is  independent  of  the  magnitude  of  a2!  Therefore, 
what  we  have  established,  is  that  all  sample  solutions  (i.e.  with  probability 
one)  (1.13)  of  the  equation  (1.11)  will  approach  zero  asymptotically.  In  the 
general  literature  on  stochastic  linear  systems,  this  is  usually  referred  to 
as  almost  sure  asymptotic  stability  (See  e.g.  [1.2]  or  [1.22]). 

Now  we  study  the  statistical  properties  of  the  solution  process  (1.13) 
in  order  to  determine  their  asymptotic  behavior.  We  concentrate  on  the  moments. 
To  illustrate  the  basic  results  here,  it  is  sufficient  to  study  the  mean  motion, 
first  moment,  E{x(t)l. 


Since  the  solution  is  known  explicitly  by  (1.13),  then  the  mean  can  be 
evaluated  as. 


-a  t 

E{x(t)  }  =  E{e 


n(s)ds 
°  >  • 


(1.16) 


The  n-process  is  gaussian,  therefore,  since  linear  operators  on  gaussian 
processes,  yield  gaussian  processes  (e.g.  see  [1.20]),  it  follows  that 

fC 

i  n(s)ds  =  N(t)  is  a  gaussian  process. 

‘  o 


We  need  merely  calculate  the  mean  aid  variance  to  obtain  the  associated 
gaussian  density  function  for  N(t),  that  is 


p(N,t) 


(N-mN(t))- 


*2-n  ^N(t) 


(1.17) 
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We  have,  by  direct  calculation 


m^t)  =  E{N(t) }  =  E{  n(s)ds}  =  E{n(s)}ds  =  0 


t  t 


o^(t)  =  E{N2(t)}  =  E{(  n(s)dsj  2}  =  E{  dx  fds  n(x)n(s)} 

\  0  /  Jo  •'o 


(1.18) 


t  t 


t  t 


•  f  f  o  —  I  - 

=  dx  ds  E{n(x)n(s)}  =  dx  ds  o  e  'T 

ft  O  *  ry  • 


=  2a2[t  +  e  _t“l] 


Since  it  is  easily  seen  by  direct  calculation,  that 


E{e^(t)}  =  exp 


then  upon  substituting  the  result  of  (1.18)  yields 


E{*(t)>.- 


(1-19) 


As  t  approaches  infinity  it  is  obvious  that  the  term  e 
determines  the  asymptotic  behavior  of  (1.19). 

2 

Therefore,  for  a  >a,  the  mean  motion  satisfies, 


(o  -a)t 


lim  E{  x  (t )  }  = 
tf00 

with  exponential  rate  of  increase.  Yet  for  any  o2  we  have  seen  that  the 
samples  approach  zero  asymptotically  with  probability  one! 


(1-20) 


-13- 


r 


i 

i 

R 

1 

S 


i 

» 

V 

v 

1 


r» 

»  • 


Again,  this  illustrates  completely  opposite  behavior  of  the  sample 
responses  and  the  mean  response  for  systems  with  random  coefficients. 

Such  non- intuitive  features  of  systems  with  uncertain  or  randomly  fluctuating 
coefficients  make  it  imperative  that  when  assuming  statistics  of  the  coefficient 
processes  the  induced  statistical  properties  of  the  response  be  carefully 
investigated  when  actual  engineering  design  considerations  are  to  be  made  for 
structures  and  their  controllers. 

Clearly,  the  statistical  measures  may  be  quite  misleading  as  to  the  nature 
of  the  true  response  characteristics.  (See  [1.23]  for  an  early  discussion  of 
this  problem,  and  [1.24]  for  the  exact  statement  relation  asymptotic  moment 
behavior  and  asymptotic  sample  behavior). 

It  follows  that  we  must  further  question  what  the  statistics  can  tell  us 
about  the  structural  response  fluctuations. 

We  will  deal  with  these  questions  along  with  the  development  of  the 
statistical  characteristics  below. 

We  will  first  look  at  the  white  noise  coefficient  case,  the  physical  noise 
coefficient  case,  and  various  approximation  schemes  in  that  order.  We  will 
concentrate  on  what  statistical  measures  (moments,  probability  densities)  can  be 
obtained  exactly,  or  approximately. 


B 

i 

B 

R 
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I . 2  Introduction  -  Time  Invariant  Structural  Equations 

This  introduction  is  concerned  with  response  prediction  in  a  time- 

invariant  system  when  it  is  not  possible  to  specify  the  system  exactly. 
Two  questions  immediately  come  to  mind.  First,  do  situations  exist 
where  it  is  m)t_  possible  to  specify  a  system  exactly  as  is  required  for 
exact  response  prediction?  Second,  does  it  make  any  great  difference  to 
the  engineer  (in  response  prediction)  if  he  does  not  have  accurate  sys¬ 
tem  specification?  It  is  instructive  to  develop  reasons  why  the  answer 
to  each  of  these  questions  is  sometimes  yes .  Two  preliminary  remarks 
are  in  order  before  proceeding. 

'..'e  know  that  in  classical  mechanics,  applied  mechanics,  structures, 
small  vibration,  in  text  books,  and  in  most  design  procedures,  the  sys¬ 
tem  is  regarded  as  known  exactly.  That  is,  masses,  stiffnesses,  and 
dissipation  are  regarded  as  known  exactly.  Response  prediction  in  these 
"ideal"  situations  is  thus  exact,  of  course.  The  possibility  that  the 
system  specification  is  in  many  practical  engineering  situations  uncer¬ 
tain  is  never  brought  up  for  consideration. 

For  any  rational  discussion  of  the  influence  of  system  uncertainty 

on  accuracy  of  response  prediction,  it  is  important  to  clearly  specify 

the  type  of  system  under  consideration  so  that  the  scope  of  the  problem 

is  kept  within  reasonable  bounds.  Vie  assume  the  system  is  holonomic 

with  generalized  coordinates  q, ,  ...,  q  ,  where  n  is  finite;  thus,  the 

i  n 

system  has  n-degrees  of  freedom.  We  next  assume  the  configuration 
q  j  =  . . .  *  q^  =•  0  is  one  of  stable  equilibrium.  We  further  assume  small 
motion  about  this  configuration  of  equilibrium,  and  take  the  dissipation 
to  be  viscous.  To  obtain  the  equations  of  motion  by  Lagrange's  method, 
we  require  the  kinetic  energy  T,  the  dissipation  function  F,  the 
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potential  energy  V,  and  the  virtual  work  6  W  of  the  external  forces.  We 
know  [  ]  that  we  may  write 


T  •  1  VVk  ’  v  ‘  I  VA 

F  =  c.,q.q.  ,  6W  =  f.(t)6q.  (1-21) 

2  jk'j  k  j 

where  sunnation  convention  is  employed  and  the  f^(t)  are  the  external 
force  components.  T  is  a  positive  definite  quadratic  form  with  constant 
rtijk  =  •  We  assume  that  F  is  positive  definite  to  ensure  only  dissi¬ 
pation  and  take  constant  c.,  =  c,  . .  Since  the  equilibrium  configuration 

jk  kj 

is  stable,  V  is  also  positive  definite  and  we  take  c.,  =  c,  (constant). 

Jk  kj 

Lagrange's  equations 


d  3T  3T  3F  3V  £  .  . 

-j—  V*"  ~  *5 —  +  nr;--  +  -r- —  =  f  .  C  t )  ,  j  =  1  ,  .  .  .  ,  n 

dt  3  q  3q.  3q.  3q  i 

2  J  j 


(1.22) 


thus  provide  the  equations  of  notion  from  which  response  prediction  are 
obtained.  Ue  can  now  discuss  how  uncertainty  enters  system  specifica¬ 
tion  in  a  significant  manner. 

First,  observe  that  in  a  deformable  structure  (system)  we  require 

an  infinity  of  coordinates  to  specify  the  configuration.  Our  choice  of 

coordinates  qj ,  . ..,  q^,  where  n  is  finite,  immediately  points  up  the 

fact  that  we  do  not  have  enough  coordinates  to  specify  the  configuration 

of  the  system  exactly  even  though  we  know  the  organization  of  the  system 

in  terms  of  members,  joints,  masses,  etc.  However,  it  is  reasonable  to 

assume  we  can  select  a  set  of  q, ,  ...,  q  that  will  serve  for  our 

1  n 

specific  purpose,  as  is  done  in  finite  element  modeling.  Thus,  we  shall 
disregard  uncertainty  in  system  specification  due  to  coordinate  choice 
in  what  follows. 

Next,  consider  V.  The  k..  (stiffness  coefficients)  will  be  calcu- 

jk 
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lated  from  geometry  ant!  material  behavior  enploying  one  of  the  standard 
methodology.  Cven  if  considerable  effort  (say,  by  finite  element  tech¬ 
niques)  is  expended  in  computation  of  these  k.^  within  the  frame  work  of 
the  selected  q^  ,  ...»  q^,  the  influence  of  joint  behavior,  for  example, 

on  the  value  of  the  k..  can  never  be  estimated  exactly.  Further, 

jk  1 

changes  in  joint  friction  due  to  corrosion  will  cause  joint  behavior  to 
slowly  change  with  time.  Since  joint  behavior,  which  determines  end- 
conditions  for  the  relevant  members,  frequently  has  a  profound  influence 
on  member  stiffness,  uncertainty  in  joint  behavior  will  produce  uncer¬ 
tain  k  .  Further,  structural  members  may  rupture  due  to  aging  or  may 
j  k 

be  partially  inoperative  due  to  assembly  and/or  manufacturing  errors  and 

remain  undetected.  N’o  matter  how  hard  we  try  to  accurately  estimate,  at 

least  some  of  the  k.,  will  actually  have  values  different  from  what  we 

jk 

estimate. 

The  entire  area  of  passive  dissipation  mechanics  is  at  best  on  an 
insecure  foundation.  Mathematical  convenience  has  dictated  the  form  of 
F  in  ( 1 . 21 )  [seeRay];to  be  specific,  the  form  assumed  produces  linear 
dissipative  terms  in  the  equations  of  notion.  Mathematical  convenience 
is  an  important  point;  however,  experimental  evidence  is  required  (and 
we  are  not  aware  of  it)  to  demonstrate  that  viscous  damping  does  produce 


physically  accurate  response  over  the  entire  frequency  range. 

Given  the  form  if  F  in  (1.21),  there  are,  in  general,  no  reliable 


techniques  for  calculating  the  c^  in  a  rational  manner.  In  the  absence 
of  large  concentrated  dissipators  (dampers),  it  is  usual  practice  either 


to  assume  the  c^  values  are  roughly  proportional  to  the  k  ^  and/or  the 

n.,  ,  where  the  proportionality  constant  is  adjusted  to  produce  dissipa- 
j  k 

tion  in  the  first  mode  of  notion  equal  to  that  observed  in  similar 


Vi' 


mmmt 
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structures,  or  to  assume  the  c  are  more  less  constant  in  value  over 

J  k 

the  structure  with  a  multiplicative  constant  adjusted  as  in  the  previous 

case.  The  c.^  arising  from  concentrated  dampers  can  be  estimated  over  a 

frequency  interval  from  test  results;  however,  accuracy  in  the  estimates 

is  usually  not  greater  than  ±50%.  In  all  events,  we  will  rarely  have 

accuracy  in  the  c^  values  even  comparable  to  that  obtained  for  the  k^« 

We  can  calculate  or  weigh  the  elements  with  as  great  an  accuracy  as 

required.  Thus,  the  physical  elements  that  enter  into  the  computation 

of  T  have  accurately  estimated  m asses.  Given  our  choice  of  the 

q  ,  ...»  q  ,  there  are  a  number  of  methods  for  calculating  the  m., 
in  j  k 

["constant  mass  matrices"].  Accuracy  in  the  values  of  the  in.,  is  obvi- 
ously  not  going  to  be  a  problem  for  a  time  invariant  system  or  for  a 
system  in  which  the  total  mass  changes  slowly  in  a  known  manner. 

Although  we  recognize  that  the  live  load  in  a  building  or  bridge  nay  not 
be  known  precisely,  we  shall  assume  on  occasion  i_n  this  report  that  the 
m  ^  jLrJL  accurately  known .  . 

V.’e  conclude  from  the  discussion  up  to  this  point  that  uncertainties 

will  exist  in  the  c.,  and  k.,  values,  with  the  former  being  larger  than 

jk  jk  o  o 

in  the  latter.  Return  to  the  first  question  - 

"do  situations  exist  where  it  is  not  possible  to 
specify  a  system  parameters  exactly  ...?" 

The  answer  is  obviously  "yes"  and  we  have  indicated  a  few  of  the  possi¬ 
bilities. 


We  now  come  to  the  second  question  - 
"does  it  make  any  difference  to  the  engineer  (in 
response  prediction)  if  he  does  not  have  accurate 
system  specification?" 

A  few  preliminary  analytical  details  will  provide  a  framework  with  which 
to  motivate  the  answer. 


Let  us  first  write  out  (1.22)  employing  (1.211  assuming  the  o  , . .  .  ,q 

1  n 


ammmmmmasSSM 


have  been  selected  so  that  the  known  (assume  the  n 


are  not  uncertain) 


kinetic  energy  has  the  form 


T 


n 

Z 

1 


*2 

qi 


q  j  +  cjkqk  *  V*  '  fj(t)  •  J  ‘  1 

Next,  introduce  the  (nxn)  matrices 


and  the  (nxl)  coltiun  vector 


n  . 


I  q  +  Kq  +  Kq  =  f  , 
where  I  is  the  (nxn)  unit  matrix. 

Let  us  now  put  the  second  order  system  in  the  first  order  form: 


x  =  Ax  +  f 


by  employing  the  substitutions 


r 

V  *  ’ 


ft 


1 


,v 

s* 

V* 


3 


a 

wVl 


!'  •>' 


where  A  is  a  (2nx2n)  matrix  and  x  and  f  are  (2nxl)  column  vectors.  We 
note  that  (1.28), when  written  out,  consists  of  2n  first  order  equations 
whereas  (1-27)  consists  of  n  second  order  equations.  The  nice  feature  of 
(1.28)  is  that  we  can  immediately  write  down  its  solution;  if  at  t  =  0, 


x  =  x  ,  then, 
o 


A(t-o) 


A( t-t ) 


x  =  e  x  +  J  er'K'~  l/f(r)dT  ,  t  >  0 


We  also  can  take  the  Laplace  transform  of  (1.28)  obtaining 

X(s)  =  {Is  -  A}_1F(s) 
where  I  is  the  (2nx2n)  unit  matrix 


(1.30) 


(1.31) 


I 


X(s)  =  L{ x ( t ) ]  ,  F(s)  =  L[ f ( t ) ] 

and  we  must  add  die  initial  condition  x  at  t  =  0.  We  can  now  indicate. 

o  ’ 

or  ploying  (1.30)  and  (1.31),  how  uncertainties  in  the  values  of  r  and  k  , 

jk  J  k 

i.e.  in  c  and  K,  can  change  the  response  significantly. 

'  e  conclude  from  (1.31)  that  the  eigenvalues  of  A  will  determine  the 
stability  of  the  systen,  since 


i  s 


{Is  -  A} 


-1  cof  { Is-A} 


. I  s-A |  (1.32) 

(where  "T"  denotes  transpose),  which  shows  that  the  eigenvalues  of  A,  as 

determined  by 


i  I 


ci  £ 


|l2s  -  A|  =  0  ,  (1.33) 

provide  poles  of  the  Laplace  transform  of  x(t).  If  we  do  not  know  the 
elements  of  the  (2nx2n)  matrix  A  precisely,  we  cannot  unequivocally 
state  the  eigenvalues  of  A  are  in  the  left-half  complex  s-plane.  Even 


-2 


purposes  of  controlling  our  original  system  orientation  we  have  the  pos¬ 


sibility  that  the  complete  system  might  at  the  worst  be  unstable  or  have 


very  poor  control  characteristics,  or  at  best  have  excellent  control 


characteristics.  The  uncertainty  of  how  the  system  would  behave  cannot 


be  tolerated  in  space  where  adjustment  and/or  repair  would  be  difficult. 


Suppose  the  eigenvalues  of  A  are  in  the  left-half  s-plane  so  we 


know  the  system  is  at  least  stable  and  suppose  the  external  forces  f 


contain  periodic  elements.  Then,  the  possibility  of  resonance 


phenomenon  comes  up.  Notice  first  uncertainties  in  the  elements  k  ,  of 

’  ’  jk 


K  means  we  are  not  certain  where  the  undamped  natural  frequencies  of  the 


system  lie;  thus,  we  cannot  be  certain  that  some  of  the  si-ple  harmonic 


components  in  the  periodic  disturbances  will  have  fre-quencies  well 


separated  from  the  undamped  natural  frequencies  of  the  system.  Second, 


uncertainties  in  the  elements  c.,  of  the  matrix  C  means  that  we  cannot 


he  certain  that  the  forced  amplitudes  of  the  system  response  will  be 


small  if  the  frequencies  of  the  simple  harmonic  components  in  the  exter¬ 


nal  forces  (disturbances)  happen  to  be  close  to  one  or  more  of  the 


undamped  natural  frequencies  of  the  system.  If  either  of  these  situa¬ 


tions  occur,  some  of  the  components  in  x  of  (1.30)  will  be  large  leading 


ultimately  to  fatigue  failures.  Examples  of  past  cases  where  these  pos¬ 


sibilities  have  occurred  and  caused  problems  readily  come  to  mind- 


We  also  have  given  specific  examples  in  the  Ceneral  Introduction 


where  uncertainties  in  parameter  values  make  a  difference  between  what 


is  predicted  to  occur,  ignoring  these  uncertainties,  and  what  actually 


occurs.  Thus,  the  answer  to  the  second  question  is  yes! 
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In  the  sections  devoted  to  time  invariant  parameters,  our  concern  is  to 
describe  how  uncertainty  in  the  parameters  influences  response  and  what  techniques 
are  currently  available  to  quantitatively  assess  this  influence.  We  shall  divide 
our  discussion  into  subsections  that  reflect  different  interests  and  methodologies. 
The  two  main  areas  of  interest  are  derived  from  problems  of  physics  and  problems  of 
engineering.  We  shall  start  with  engineering. 

The  most  important  topic  for  engineering  systems  is  how  uncertain  parameter 
values  influence  the  accuracy  of  system  response  prediction.  It  often  suffices 
to  know  how  these  uncertainties  influence  the  accuracy  in  estimating  the  values 
of  the  natural  frequencies  and  their  corresponding  normal  modes  of  motion  in  a 
conservative  system  (C=0).  Since  linear  systems  response  prediction  depends 
upon  frequency  response  or  impulsive  admittance  (i.e.  Green's  function,  impulse 
function),  our  interest  centers  on: 

natural  frequencies 
normal  modes 
frequency  response 
impulse  response 

In  some  situations  we  will  be  directly  interested  in  the  response  q. 

It  is  important  to  note  that  there  are  two  ways  to  quantitatively  characterize 
uncertainty  in  the  parameters.  We  may  simply  have  bounds  on  the  parameter  values. 
Alternatively,  we  may  consider  the  parameters  as  random  variables,  described  by 
their  joint  probability  distributions,  or,  at  least,  by  their  first  two  moments. 

We  shall  treat  both  of  these  characterizations,  since  each  can  arise  in  application. 
The  broad  classes  of  techniques  available  to  pursue  these  subjects  are: 

perturbation  methods 
Liouville's  equation 
mean-square  approximate  systems, 
bound  determination. 

We  shall  discuss  each  of  these  techniques  separately,  realizing  that  there 
will  be  overlap. 


II.  Moments  -  White  Noise  Coefficient  Case 

In  this  section  we  shall  be  concerned  with  the  case  that  the  coefficients 

of  linear  differential  equations,  contain  white  noise  components.  In  particular, 

we  are  concerned  with  Gaussian  white  noise  coefficients.  The  white  noise 

{W  ,  tCHO,00]},  is  characterized  by  the  fact  that  its  power  spectral  density  is 

constant  over  the  entire  frequency  domain  (-“,“).  Therefore,  by  the  Fourier 

transform  relation  between  power  spectral  densities,  f(w),  and  covariance 

functions,  y(x),  it  follows  that  the  covariance  for  the  white  noise  is  an  impulse 

function.  Thus  for  the  white  noise,  f  (w)  5  S  ,  and  its  Fourier  transform 

w  o 

(the  covariance)  is  y(x)  =  2ttS  6(x).  It  immediately  follows  from  y  (x),  that 

o  w 

W(t),W(t+x)  are  uncorrelated  random  variables,  no  matter  how  small  x.  Further- 

2 

more  since  formally,  we  must  have  E{w  (t)} =  2xtS  6(0)  which  is  undefined,  we  see 

o 

that  the  white  noise  would  possess  infinite  power,  making  it  a  mathematical 
concept  rather  than  the  model  of  a  process  that  is  observed  in  nature.  Moreover, 
if  we  assume  that  W(t)  is  a  Gaussian  random  variable,  then  it  follows  that  W(t), 
W(t+x)  are  not  only  uncorrelated,  they  are,  indeed,  independent  random  variables. 

The  question  that  comes  up  is,  how  do  we  interpret  the  meaning  of  differential 
equation  models  of  real  systems,  with  Gaussian  white  noise  coef f icients?  It  was 
the  classic  work  of  K.  Ito  [2.1],  in  the  1940's  that  answered  this  question,  long 
before  it  was  of  importance  to  modern  optimization  and  control  applications. 

Although  a  thorough  development  of  these  ideas  is  beyond  the  scope  of  this 
report,  we  shall  illustrate  a  few  of  the  basic  ideas.  The  interested  reader 
should  consult  texts  such  as  [2.2],  [2.3]  for  the  fundamental  development.  The 
key  point  is  that  the  Gaussian  white  noise,  W(t),  possesses  a  representation  via 
the  Brownian  motion  process  (B(t),  t£[0,°°]}. 

An  especially  readable  account  of  the  interpretation  can  be  found  in  [2.4] 

Vol .  II  ,  as  well  as  the  recent  book  on  oarametric  excitation  bv  Ibrahim  [2.51 . 

The  Brownian  Motion  is  discussed  in  almost  all  elementary  tests 
on  stochastic  processes. 

For  0<t1  <t2»--<tn,  for  any  {t^}  and  n ,  the  joint  density  is 
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The  Brownian  Motion  process,  also  referred  to  as  the  Wiener  process  after 
Norbert  Wiener  who  first  made  a  thorough  study  of  the  properties  in  1923  [2.6], 
and  1930  [2.7],  is  defined  by  the  class  of  Gaussian  joint  probability  densities 
(2.1).  The  Brownian  Motion  process  is  a  Gaussian  process  with  stationary, 
independent  increments  satisfying. 


P{B(0) =  0} =  1,  E{B(t) } =  0,  E{B(s)B(t)} =  y  ( s ,  t ) =  a 2min ( s , t )  (2.2) 

D 

Among  Wiener's  fundamental  contributions  to  the  development  of  this  process, 
was  to  determine  that  the  sample  functions  of  the  process  are  continuous 
functions  on  any  finite  interval,  and  that  the  sample  functions  are  nowhere 
differentiable. 

Now  from  (2.2),  we  see  from  the  covariance,  that  the  Brownian  Motion  is 
non-stat ionary.  Furthermore,  from  elementary  properties  of  derivatives  of 
processes,  it  follows  that  the  covariance  of  B(t),  is  given  as 

-2  , 

J  YB(s,t) 

os3t 

This  result  is  immediately  obtained  formally,  from 


l2Vs’t:> 

^sSt 


7 

)x")t 


Et  B(s)B(t)  / 


=  El B(s)B(t)  }  . 
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Upon  taking  these  derivatives  of  yD(s,t)  given  in  (2.2),  we  also  obtain  the 

O 

formal  relation, 


3  YR(s,t) 

-JZTt — =  0  5(t"s)  =  Yw(t_s)  • 


(2.3) 


Thus,  we  obtain  the  representation  for  the  Gaussian  white  noise  as. 


B (t)  -  W(t). 


(2.4) 


Therefore,  for  systems  with  Gaussian  white  noise  coefficients,  we  can  replace 
the  white  noise  terms  with  the  formally  differentiated  Browmian  motions. 

We  write  the  general  linear  differential  equation  with  white  noise 
coefficients  as, 


x(t)  =Gx(t)  +  I  HHj  (t)  x  (t)  . 
i  =1  1 


(2.5) 


where  the  nxn  matrices  G.CH^}  are  known,  and  the  (W^(t)]^_^  are  independent 
white  noise  coefficients,  and  x  is  an  n-vector.  However,  since  we  can  represent 
W^(t)  =B^(t),  then  (2.5)  may  be  rewritten  as, 


x(t)=Gx(t)+[  H1B.(t)x(t)  . 

i=l  1 


(2.6) 


As  we  stated  above,  B^(t)  is  only  a  formal  derivative,  since  as  a  result  of 
Wiener's  investigations,  the  Browmian  motion  does  not  possess  derivatives.  The 
approach  of  K.  Ito,  was  to  interpret  (2.6)  as  an  equation  in  differentials, 


K 

dx(t)  =  Gx(t)dt  +  j  H1dB.(t)x(t)  . 


(2.7) 


The  meaning  of  this  equation  is  via  the  integral  equation. 


x(t)-x(a)=  Gx(s)ds  +  7  H  1x(s)dB .  (s) , 

•'a  i=l  U  1 


(2.8) 
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where  the  Integral  in  (2.8)  is  the  so-called  Ito  stochastic  Integral  and  the 
equation  (2.7)  is  generally  referred  to  in  the  literature  as  the  Ito  differential 
equation. 

Thus,  upon  starting  with  (2.5)  as  the  formal  differential  equation 
describing  a  system  with  white  noise  excitations,  we  are  lead  to  the  math¬ 
ematically  meaningful  form  (2.7). 

The  importance  of  this  form  lies  in  the  many  properties  of  the  Brownian 
motion  differential  coefficients  dB(t).  In  particular,  since  the  Brownian 
motion  process  is  a  process  with  independent  increments,  then  for  each  i, 
dBi(s),  dBi(t)  are  independent  random  variables  for  s  ^  t.  (We  define  the 
differential,  dB(t) =  B(t +  dt)-B(t)  for  positive  time  differentials  dt>0). 

Also,  as  a  result  of  this  interpretation  of  the  differential,  it  also  follows, 
that  dB^(t)  is  independent  of  x(t)  for  every  i.  We  can  interpret  this 
formally  as  follows.  Through  the  differential  equality  (2.7),  x(t)  is  a 
functional  of  (B^(s),  sit)  for  all  i.  However,  since  the  Brovnian  motion  is 
an  independent  increment  process,  it  follows  that  dB^(t) = B^(t +  dt)-B^(t)  is 
independent  of  all  combinations  (or  functionals)  of  (B^s)  sit).  This 
independence  is  of  fundamental  significance.  We  can  see  the  immediate  effect 
of  this  property,  when  we  investigate  the  mean  value  E{x(t)},  for  the  solution 
response  to  (2.7). 


By  taking  expectations  directly  on  (2.7)  we  obtain. 


dE{ x(t) } =  G  E(x (t) }dt  +  l  R1E{dB.(t)x(t)}. 

i=l 


(2.9) 


Due  to  the  independence  of  x(t)  and  the  dB^(t),  we  can  write  using  the  fact 
that  EfdB(t)} =  0. 


E{dBi(t)x(t) } =  E{dBi(t)}  E{x(t)} =  0 


(2.10) 


Therefore,  we  immediately  see  that  the  mean  equation  from  (2.9),  (2.10)  is 


—  E{x (t ) }  =  G  E{x(t)  } , 


(2.11) 


which  is  immediately  solvable  as  a  linear  vector  equation  with  constant 
coef f icients. 
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Because  of  these  properties,  all  moment  statistics  of  the  solution 
process  to  (2.7)  may  be  obtained.  Perhaps,  the  most  important  result  of  the 
properties  of  the  Brownian  motion,  and  the  form  of  the  Ito  differential 
equation,  is  that  the  solution  process  [x(t),  tglt^,®)},  is  a  vector  Harkov 
process.  See  [2.3]  and  [2.12]  for  discussion  of  this  deeply  significant  fact. 

Moreover,  it  was  the  motivation  to  represent  a  Markov  process  explicitly 
that  lead  K.  Ito  to  study  the  differential  equations  that  bear  his  name. 

The  important  tool  by  which  we  can  set  up  the  various  moment  equations 
is  the  so-called  generator  or  backward  operator  for  the  Markov  solution 
process  of  (2.7).  We  will  obtain  this  generator  via  the  characteristic 
functional,  an  approach  that  appeared  in  an  important  early  article  by 
Moyal  [2.8]  and  is  attributed  to  the  statistician  M.S.  Bartlett.  This 
approach  will  be  used  to  obtain  the  Liouville  equation  in  our  review  of  the 
random  constant  coefficient  case,  following  the  approach  in  [2.91.  These  ideas  were 
also  applied  in  a  fundamental  paper  on  statistical  turbulence  theory  by 
E.  Hopf  [2.10]  and  later  by  the  physicist  S.F.  Edwards  2.1 

To  illustrate  these  ideas,  we  write  the  development  for  the  scalar  case. 

The  n-th  order  operators  can  then  be  written  down  immediately.  For  the  simple, 
scalar  Ito  equation  (general  linear  or  non-linear)  we  write. 


dxt  =  g(xt)  dt  +  h(xt)  dBt. 


(2.12) 


We  are  interested  in  the  characterstic  functional 


Jx(t,u).E{eto<t)) 


(2.13) 


Upon  taking  the  differential  with  respect  to  t  of  (2.13),  we  will  obtain,  upon 
a  formal  interchange  of  expectation  and  differential  operators  - 


dt?x(t’u) 


E(dteiux(t)> 


=  E{[iudx(t)  +  y  (iua  x(t) )2  +  o(dt) ]  eiuX<'t^} 

=  iuE{dx(t)  eiux(t)]  +|  (iG)2  Ef(dx(t))2  eiux(t)} 


(2.14) 


+  o (dt) 


WWW 
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Here,  we  have  taken  the  first  two  terms  in  the  expansion  since,  as  we  shall  see 


from  (2.13)  there  will  be  a  contribution  from  the  second  order  term.  Higher 


order  terms  (dx(t))  ,  k>2  will  all  be  o(At),  and  thus  will  not  yield  a 


contribut ion. 


From  (2.13),  we  have 


d^x(t,u)  *  iu  E{  [g(x(t)dt  +  h(x(t))dB(t)  ]  eiux^t^ } 


(2.15) 


+  |  (iu)2  E{[(g(x(t))dt)2  +  2g(x(t))h(x(t)dtdB(t) 


+  h2(x(t))  (dB(t))2]  eiux(t) } 


+  o  (dt) 


Now,  since  dB(t)  is  independent  of  any  function  of  x(t),  as  discussed  above, 
we  have  from  (2.2),  for  any  f(x(t)). 


E{dB (t)  f(x(t))}  =  E{dB(t) }  E{f(x(t)}  *  0 


(2.16) 


E{(dB(t))2  f (x(t) }  -  E((dB(t))2}  E{ f (x (t) } =  a  dtE{f(x(t)} 


Therefore,  we  can  write  (2.15)  as, 


dt»x(t,u)  =  iu  E(g(x(t))  eiux(t)}  dt 


+  ~  (iu)2  E{h2 (x (t) )  eiUx(°}  dt 


(2.17) 


+  o (dt) 


Finally,  we  obtain  upon  dividing  by  dt,  and  letting  dtl-0, 


51?  (t,u)  n2  2  2  iux(t). 

-  -E{[iug(x(t))+~  (iu)2  h2(x(t))]  e  UxU)} 

7  t  A 


(2.18) 


ryt 


We  must  now  recall  the  elementary  fact  that 


*fx(t,u)  =  |  eiu  p(x,t)  dx, 

—  QO 

that  is  the  characteristic  function,  and  the  corresponding  probability  density 
function  are  Fourier  transform  pairs. 

Upon  taking  inverse  transforms  of  the  equation  (2.18),  using  the  operator 
identity  (iu)  -  — j-  ,  (2.18)  yields  the  partial  differential  equation  for 

JX 

p(x.t) , 


■> P ( x ,  t )  _  _  3  [g(x)p(x,t)  ]  0^  32  [h2(x)p(x,t)  ] 

3t  3x  2  2 

3x 


(2.19) 


We  recognize  this  equation  hs  the  Fokker-Planck  equation  that  has  been  so 
important  in  the  study  of  Markov  diffusion  processes. 


The  operator 


-3[g(x)  •  ]  of  32[h2(x)  ♦  ] 
3X  2  3x2 


(2.20) 


is  referred  to  as  the  Forward  operator  in  the  literature.  See  [2.12]  as  well 
as  other  sources. 

The  adjoint  operator 


5  8(x)  ^  +  Y  h2(x)  ^2  •  (2.21) 

3x 

is  referred  to  as  the  generator  of  the  process  defined  by  the  Ito  differential 
equation  (2.12). 

The  generator  is  the  most  important  operator  for  diffusion  processes 
since  it  always  exists.  It  is  known  that  the  forward  operator  may  not  exist 
for  certain  Markov  processes.  This  is  an  advanced  concept  which  is  usually 
discussed  in  fundamental  studies  of  diffusion  processes  [See  e.g.  [2.13]]. 


-28- 


It  is  the  generator  that  is  of  interest  to  us.  It  is  easily  seen  that 
for  any  f(x),  from  (2.19),  since  we  may  write. 


^  E{f(x(t))}  -  ~  J  f(x)p(x,t)  dx  -  f(x)  dx. 


we  obtain  the  equation  for  the  mean  value  of  f(x(t), 


±  E{f  (x  (t)  )  }  =  E^f(x(t))} 


(2.22) 


The  equality  (2.22)  is  a  simplified  form  of  what  is  referred  to  in  the 
literature  as  Dvnkin's  theorem,  [see  [2.13],  [2.14]].  The  apparent 

first  application  of  these  ideas  to  structural  systems  was  presented  in  [2.15], 
where  stability  of  the  second  moments  of  the  linear  oscillator  with  white  noise 
coefficients  was  studied. 

k 

Considering  only  the  moments  ul  (t)  *  E{x  (t)},  we  have  from  (2.21),  (2.22), 


~  mk(t)  =  kE{g(x(t)xk  E(h2(x(t))xk‘2(t)} 


(2.23) 


It  is  clear  that  this  equation  cannot  be  solved  for  arbitrary  functions 
g(x),h(x).  But,  to  our  good  fortune  the  linear  case  is  completely  determined 
by  (2.23).  Thus,  for  the  linear  form  of  (2.12)  where  g(x)  *  ax,  h(x)=*bx,  the 
moment  equation  (2.23)  becomes 


m^t)  -  kE'rax(t)x(k~1)(t)}^  —  E(b2x2(t)x  (k_2)  (t) } 


which  the  reader  can  easily  put  into  the  form, 


d  r,  ,  k(k-l)  2,2, 

—  ra^ft)  =  [ka  +  — - -  3  b  ]  riL(t), 


(2.24) 


which  yields. 


E  (t )  =  m^(t )  =  m  (0)e 


for  the  k  moments. 


In  fact  for  any  n  order  linear  Ito  equation  we  obtain  a  consistent  set 
of  equations  for  the  ktb  moments  that  can  be  solved  exactly.  We  see  this  by 
first  presenting  the  generator  (backward  operator),  for  the  general  linear 
ntb  order  Ito  equation  (2.  7  ).  It  is  obtained  directly  as 

(ISiiV  ^  bij(X)  3x\x  1 

x  i-1  13  3  j  i,j-l  3  i  j 


where  by  (x)  -  ^  a\ 


a*  dt  »  E{  (dB£(t))2h  and  -  (h^  )  . 

We  note  that,  exactly  as  for  the  first  order,  scalar,  example  above, 

the  coefficient  of  the  -r —  term  is  linear  in  x,  and  the  coefficient  of  the 

3Xi 


3x .  3x , 


term  is  quadratic  in  x. 


Hence,  there  will  always  be  a  closed  set  of  equations  to  solve  for  the 


i  th 

k  moments. 


In  fact,  we  immediately  have  for  the  expectation  of  the  general  function 
f(x^(t),... ,xn(t)), 

—  E{f  (x,  (t)  , . .  •  ,x  (t))>  *  E{«S?f  (x.  (t) ,. . .  ,x  (t))}  *  (- 

dt  1  n 

similarly  to  (2.22). 

We  look  at  the  classic  example,  apparently,  first  studied  in  [2.15]. 

This  Is  the  second  order  linear  oscillator,  with  a  white  noise  coefficient. 


x(t)+2ijLx(t)  +  <w  +W(t))x(t)  -  0  (2 

Putting  <2.21 ,  into  a  linear  Ito  form,  we  set  =  x  ,  X2*x  and  noting 
that  dB(t)  *  W(t)dt ,  as  in  (2.  7),  the  second  order  oscillator  equation  becomes 


dx^  “  x0dt 


'dx_  •  -(2'.x0  +  ~“x  )dt  -  :<  dB(t) 

L  k.  L  i 


mmm 


H 


The  generator  of  the  (x^jX^)  diffusion  process  can  immediately  be  written 


i? 


3  .  2  3  a2  2  32 

X2  3^  -<2<“x2  +  W  Xl)  ^2+“  X1  ^ 


(2.29) 


2 

For  the  three  moments  m2^(t) =  E{x^(t)>,  m^2 (t ) *  E{ x^ (t )x2 (t)  , 
m22 (t)  =  E^x2 (c) ^ »  the  relation  (2.26),  for  J^given  by  (2.29)  yields  the  equations. 


/rail(t) 

31  I  "l2(t) 


m22 


-2u  -4nw/  \n»22(t) 


»u(t) 

"l2(t) 


(2.30) 


which  can  be  solved  explicitly. 


i  I 

I, 


For  the  general  n  order  system  (2.7),  one  can  easily  show  that  the  second 

moments  E{x  (t)x  (t) }  =  m  (t),  can  be  expressed  as  the  solutions  of  the 
u  v  uv 

differential  equations 


K  n  „  „ 

l  az0  l  h  h  m  (t) 

a-i  *  r ,s,=i  ur  vs  rs 


(2.31) 


where  u,v=l,2,...,n. 

Here,  we  have  used  (2.25)  directly. 

Finally,  we  shall  mention  a  generalization  of  the  2nc*  moment  formula 
(2.31)  that  holds  for  all  pfc^  order  moments,  [2.16],  [2.17]  for  the  linear 
Ito  equation  (2.7). 


Motivated  by  the  algebraic  theory  of  linear  differential  equations, 

one  can  define  for  a  given  n-vector  x  and  a  given  positive  integer  p,  the 

[p] 

associated  vector  x  whose  components  are 


/.  -■  «-  /  .*  .*  .* 

.s 
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The  vector  satisfies  ||  x^  ||  =  ||  x||P  ,  where  j|x  =  (x,x)  and 

more  generally,  (x,y)P  =  (x^,y^).  Furthermore,  this  concept  extends  to 
matrices  through  the  definition  y  ■  Ax-  One  defines  ,  as  that  matrix  that 
yields 


y[p]=A[p]x[p] 


(2.34) 


For  linear  systems 


x  (t ) =  A(t)x(t) 


(2-35) 


the  differential  equality 


x(t+h)  =  (I  +hA(t))x(t)  +  0(h) 


holds,  yielding  from  the  definition  (2.34) 


:Ip](t+h)  =  (I  +  hA(t))[p]xtp](t)  +  0(h2) 


(2.36) 


Upon  defining  the  limit 


lim  i  ((I  +  hA(t))[p]-  I[p1)) 
h-t-o 


k[p]  ’ 


we  now  determine  the  associated  differential  equation  for  , 


i[p] (t)  =  A(t)  x[p] (t) 


(2.37) 


In  order  to  relate  these  ideas  back  to  the  original  Ito  equation  (2.7), 
one  can  apply  the  usual  Ito  calculus  [2.12]  to  obtain  the  Ito  differential 
equation  for  x^  as 


aasamaaiii^^ 
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Applying  expectations  to  the  Equation  (2.38)  immediately  yields  the 
linear  differential  equation  for  the  pth  order  moments  as 


[p] 

=  A  E{x[p](t)}  .  (2.39) 

P 


The  result  (2.39),  can  also  be  obtained  directly  from  application  of 
(2.2')  f or  <£  given  by  (2.25). 

Although,  in  principle,  the  equations  for  the  pth  moments  for  linear 
homogeneous  Ito  differential  equations  have  been  known  for  many  years,  the 
general  form  (2.39)  can  be  quite  useful. 

These  ideas  extend  quite  simply  to  the  non-homogeneous  case  as  well. 

We  mention  that  the  results  (2. 32)-(2.39) ,  have  been  derived  from  algebraic 
considerations  and  more  particularly  by  consideration  of  the  Lie  algebras 
[see  [2.16)]  generated  by  the  matrices  (G,  H^, 

Finally,  it  should  be  noted  that  there  is  no  conceptual  or  analytic 
difficulty  to  obtain  the  moments  for  the  linear  system  with  external  forcing 
functions.  The  general  linear  Ito  equation  (2.7)  would  become, 

K 

dx(t)  =  Gx(t)dt  +  l  H1dB.(t)x(t)  +  FdV(t), 
i=l  1 


where  F  is  a  constant  nxm  matrix,  the  vector  V  is  an  m-vector  of  Brownian 
motions  that  are  independent  of  the  (B^(t)},  and  are  independent  among  their 
components. 


I 


,m*j 


■  o 
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For  this  case,  the  generator  for  (2.40)  would  be  +EEV>  where  EE  is 

X  A  X 


given  by  (2.25)  and 


CP  1  r  3^ 

*r  =  2  ^ .  aij  3x3x.  * 

i.J=l  1 


where  £„£.,  , 


and  F  =  ( f  ) 
nxm  ij 


The  equation  (2.40)  is  a  non-homogeneous  equation.  The  associated  moment 

equations  that  are  derived  from  the  addition  of  EE  given  by  (2.41)  to  EE  given 

F  x 

by  (2.25)  would  also  generate  a  set  of  non-homogeneous  linear  equations.  Finally, 
we  wish  to  comment  on  the  probability  densities  for  (2.7).  We  know  that  the 
probability  densities  for  the  solutions  to  Ito  differential  equations  are  the 
solutions  to  the  Fokker-Planck  equations,  given  by 


’‘EE  p,  where 


(2.42) 


EE  is  the  adjoint  of  the  generator 

Thus,  for  example,  for  the  general  linear  equation,  the  adjoint  otE£^ 
given  by  (2.25)  is 


>8&-- 

x  ,L ,  3x.  2  ,L .  ,  3x3x. 

1=1  i  i,J=l  J 


(2.43) 


For  the  simile  scalar  equation 


dx(t)  =  Bx(t)dt  +  xdB(t),  E{(dB(t))^}  =  a^dt. 


(2.44) 


msHsmm 


w 
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whose  Fokker-Planck  equation  is 


3p(x,t)  _  a  3 [xp (x »t) ]  of  a2[x2p(x,t 
3t  3x  ^  2  ,  2 

JX 

we  can  obtain  the  probability  density  explicitly  as. 


[log  x  -(6  -  a2)t]2/2a2t 


P(x,t)  = 


(2.45) 


Unfortunately,  for  higher  order  linear  systems,  one  cannot,  in  general, 

solve  for  the  probability  densities.  In  certain  cases,  however,  the  stationary 

density  as  the  solution  of 


cp*  / 

=  p(x) 


(2.46) 


can  be  obtained. 
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III.  Approximately  White  Noise  Coefficients 

As  we  discussed  in  II,  a  white  noise  is  distinguished  by  the  fact 
that  its  spectral  density  is  constant  over  the  entire  frequency  domain  (-<*>,  <*>) . 

Although  the  white  noise  itself  is  a  mathematical  abstraction,  not 
present  in  nature,  there  are  real  or  "physical"  noise  processes  that  are  wide 
band.  These  processes  may  possess  power  spectral  densities  that  are  essentially 
flat  over  a  broad  frequency  range  and  then  exhibit  a  rapid  drop-off  to  neglible 
frequency  content. 

It  has  been  a  common  procedure  throughout  the  development  of  stochastic 
methods  to  deal  with  problems  of  random  excitations  to  replace  such  wide  band 
processes  with  White  noise. 

In  the  case  that  the  wide  band  gaussian  physical  noise,  n(t) ,  is  an 
external  excitation,  such  as  in  the  simple  oscillator, 

x(t)  +  2^ujx(t)  +  (o2x(t)  =  n(t) ,  (3.1) 

the  typical  procedure,  over  the  years,  has  been  to  replace  n(t)  by  W(t) 

(gaussian  white  noise)  and  proceed  with  the  analysis  to  obtain  the  solution 
process,  moments,  probability  densities,  etc.  (In  this  linear  case,  (x,  x) 
will  be  gaussian  random  variables  for  gaussian  excitations). 

For  external  excitations  this  procedure  can  be  justified.  Difficulties 
occured  when  researchers  in  random  vibrations  attempted  to  make  the  same  type 
of  replacements  for  wide  band  random  coefficients.  In  the  random  coefficient 
case,  the  replacement  cannot  be  simply  made.  A  deeper  analysis  is  required. 

For  the  early  discussions  and  ultimate  clarification  of  these  questions,  the 
reader  is  referred  to  [3.1],  [3.2],  [3.3],  [3.4],  [3.5]. 

The  basic  problem  can  be  illustrated  via  the  first  order  differential 
equation. 


x(t)  +  n(t)  x(t)  =  0. 


(3.2) 


I 


r  ’  «k  *  .*  >  O* 
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For  any  "physical"  noise  n(t),  whose  sample  functions  are  well  behaved 
in  the  sense  that  they  are  at  the  very  least  piecewise  continuous,  the  solution 
may  be  represented  as, 

-[  n(s)ds 


x(t)  =  x  e  0 
o 


(3.3) 


where  x(0)  =  x  . 

o 

The  natural  question  that  occurs  is  as  follows: 


"If  n(t)  is  a  wide  band  gaussian  process  and  we  replace  n(t)  by 
the  gaussian  white  noise  W(t)  in  (3.2),  will  the  solution  be  given  by  (3.3) 
where  n(s)  is  replaced  by  W(s)  in  the  integral?" 


This  is  exactly  what  was  done  in  the  early  1960's  in  order  to  study 
oscillators  with  white  noise  coefficients.  This  would  allow  us  to  write  (3.3) 
as 


x(t)  =  x  e 
o 


-j  W(s)ds 

'  r\ 


dB  (s) 


=  x  e 
o 


=  x  e 
o 


-B(t) 


(3.4) 


from  the  representation  (2.4)  of  the  gaussian  white  noise  in  terms  of  the 
Browmian  motion.  We  would  have  to  verify  that  the  sample  solution  (3.4)  satisfies 
the  original  equation 


x(t)  +  W(t)  x(t)  =  0, 
or  in  the  proper  Ito  differential  form 

dx(t)  +  x(t)dB(t)  =  0  (3.5) 

(again  identifying  W(t)dt  -  dB(t)). 

This  is  simply  obtained  by  taking  the  differential  of  the  solution  in  (.3.4), 

-B(t) 

x  e 
o 

But,  differentials  of  functions  of  Browmian  motion  possess  a  somewhat 
different  form  than  differentials  of  the  ordinary  calculus. 


s 
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This  was,  in  fact,  a  major  point  in  the  development  of  Ito's  stochastic 
differential  equations  [see  eq.  [3.6],  [3.7]].  Essentially,  since 

2  o 

E{(B(t  +  At)  -B(t)  }=  a*At,  for  any  At>o,  it  is  known,  on  a  sample  property  basis 

2 

that  (B(t  +  At)  -B(t))  =  a2At  with  probability  one.  Therefore,  upon  taking 

a  Taylor  expansion 

2 

dF  (B  (t) )  =  F'  (B(t)dB(t)  +  F"(B(t)) 


+  ...  +F(n>(B(t))-^p-n+  ... 


(3 . 6) 


2 

we  must  keep  the  (dB(t))  ~  a2dt  term,  all  higher  order  terms  (dB(t))  -  o(dt) 

for  k>2. 

This  leads  to  the  differential  formula, 


dF(B(t) )  =  F'(B(t))dB(t)  +  j-  F"(B(t))dt, 


(3.7) 


which  is  referred  to  in  the  literature  as  the  Ito  differential  formula. 

This  is  the  proper  relation  that  must  be  used  to  study  the  calculus  of  stochastic 
differential  equations  with  gaussian  white  noise  coefficients.  For  more  general 
functions  F(B(t),t),  one  can  show  that 

o2 

dF(B(t),t)  =  Ftdt  +  FBdB+|  FRBdt  .  (3.7a) 

Upon  application  of  the  formula  (3.7)  to  the  function  (3.4),  we  find  that 
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Thus,  simply  by  replacing  n(s)  with  W(s)  in  (3.3)  does  not  yield  the  proper 
solution  to  (3.2)  with  n(t)  replaced  by  W(t). 

n(s)ds 

Thus,  XQe  *  o  is  the  solution  to  an  equation  that  differs  from  the 

ft 

W(s)ds 

equation  for  which  xQe  '°  is  a  solution.  Yet,  for  the  ordinary  deter¬ 

ministic  calculus,  the  equation  would  be  the  same.  Thus,  we  must  look  further 
at  how  to  study  the  wide  band  coefficient  case  by  replacement  with  white  noise 
coefficients.  It  appears  from  our  simple  example  that  the  equation  would  have 
to  be  modified  to  yield  the  same  analytical  results. 

We  can  show,  simply,  how  the  equation  must  be  modified.  Consider  the 
scalar  equation, 


dx(t)  =  f (x(t) )dt  +  g(x(t))  dy(t) ,  (3.9) 

where  C ^  exists  as  a  physical  noise  process. 

We  assume  that  we  can  write  the  solution  as 


x ( t )  =  F(y(t),t)  (3.10) 

Thus,  we  would  find 


dx(t)  =  F  (y(t),t)dt  +  F  (y(t) ,t)dy(t)  (3.11) 

y 

Therefore,  we  must  have  by  identifying  terms  in  (3.9),  (3.11), 

Ft(y(t),t)  =  f (F(y (t ) , t) ) 

F  (y(t),t)  =  g(F(y (t) ,t) ) .  (3.12) 

Now  suppose  we  replace  y(t)  by  B(t)  in  (3. 12), we  set  a  new  x(t)  equal  to, 

x ( t )  =  F(B ( t ) , t )  .  (3.13) 
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Upon  applying  the  Ito  differential  formula,  (3.7a)  we  will  have 


dx(t)  =  Ft(B(t),t)dt  +  FB(B(t),t)dB  +  \  FBB(B(t),t)dt 


(3.14) 


now  F_(B(t),t)  =  g(F(B  (t) , t) ) .from  (3.12).  Therefore, 

D 


FBg(B(t),t)  =  g ' (F(B (t ) , t) )  FB(B(t),t) 


-  g'(F(B(t),t))  g(F(B(t),t)), 


(3.15) 


again  from  (3.12). 


Upon  realizing  that  (3.13)  holds,  then  the  substitutions  of  (3.12),  (3.15) 
into  (3.14)  will  give  us  the  equation 


dx(t)  =  f (x(t)dt  +  g(x(t))dB(t)  +  —  g' (x(t))g(x(t)dt 


(3.16) 

We  see,  therefore,  the  change  in  the  equations  (3.9),  (3.16)  through  the 


term 


g'(x(t))  g(x(t))dt. 


(3.17) 


Therefore,  it  is  just  this  term  that  must  be  added  to  the  equation  (3.9) 


when  replacing  the  coefficient  by  the  Gaussian  white  noise  W(t)  | 


The  term  (3.17)  is  a  correction  term  and  is  usually  referred  to  in  the 
literature  as  the  Wong-Zakai  correction  term. 


For  the  general  n-dimensional  vector  differential  equations,  we  will  find 


that 


dx(t)  =  f(x(t),t)  dt  +  G(x(t) ,t)dy(t) , 


(3.18) 


where  x,f  are  n-vectors,  y  is  an  m-vector  and  G  is  an  nxm  matrix,  will  be 
replaced  by  the  equations 
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dx±(t)  = 


[2  n  m 
f  (x,(t),t)  +  —  l  l 

k=l  1=1 


3g  (x(t) ,C) 

gk£(x(t)tt)  - ^ -  dt 


+  [  g . „ (x(t),t)dB  (t)  ,  i=l , . . . , n 


(3.19) 


The  immediate  question  for  us,  is  how  will  this  change  the  linear 
equations  with  random  coefficients  that  are  "almost"  white  noise.  We  have 
already  seen  how  the  simple  first  order  equation  (3.2)  is  modified  to  obtain  the 
added  term  present  in  (3.8).  Our  interest,  of  course,  is  for  higher  order  systems. 

Thus,  we  can  for  example  consider  the  case  of  the  second  order  oscillator 


which  corresponds  to 


dx  =  x  dt 

^  dx^  =  -(2£idX2  +  w  x^)dt  -  x^c 


(3.20) 


x-f  2£cux  +  (uj  +y)x  =  0, 


where  y  is  a  physical  noise  coefficient. 

Identifying  terms  in  (3.18),  (3.19)  with  the  system  (3.20),  we  see  that 


G(x)  = 


(3.21) 


Hence,  the  correction  terms  in  (3.19)  become 


Z  .  ^ 

A  Skl  3x. 
k=i  k 


il  _ 


i  =  2 


(3.22) 


mm® 
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Thus,  there  is  no  correction  required.  The  associated  Ito  equation  is 

2 

x  +  2so)X  +  (u  +B)x=  0 

Indeed  for  the  second  order  oscillator,  a  correction  term  is  required 
only  when  there  is  a  randomly  fluctuating  parameter  in  the  damping  term. 

Indeed,  for  the  n*"^1  order  general  linear  system,  there  will  be  a  correction 

t  h 

term  required  only  when  there  is  a  random  parameter  appearing  in  the  (n-1) 
derivative  te'ms.  No  corrections  will  be  required  to  the  linear  system  equations 
when  the  randomly  fluctuating  coefficients  appear  in  terms  lower  than  the  (n-lj"*1 
derivative.  The  following  example  illustrates  this  property. 

In  exactly  the  same  fashion  as  above,  we  can  easily  see  that  for  the 


dx1  =  x^dt 


dx2  =  ~(2r,bix^  + d  x^)dt  -  x^dy 


(3.23) 


corresponding  to 


x  +  (2cw+y)x  +  a)  x  =  0, 


the  G  matrix  is  given  as 


’21  /  . 


(3.24) 


Hence,  the  correction  terms  become. 


, 1  .  8kl  •x, 
k=  L  k 


(3.25) 


Thus,  the  associated  Ito  equation  becomes 


dv  =  x?dt 


dx^  =  [  —  ( 2 '  +  j‘x^)  +'V'  xo  ]dt  -  x^dB, 


(3.26) 
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which  corresponds  to 


x-t-  (2',^  -  ^"  +  B]x  +  ^  x  =  0. 


(3.27) 


Notice,  that  for  the  associated  Ito  equation  the  -  —  x  has  a  destabilizing 
effect  on  the  overall  system. 

Naturally,  we  easily  see  how  these  correction  terms  will  appear  for  the 
linear  system.  The  important  application  is  that  of  determining  the  approximate 
moments  for  the  specific  system. 

For  the  physical  noise  coefficient  system  (3.20),  we  cannot  determine  the 
various  moments. 

However,  for  the  associated  Ito  system,  with  no  correction  terms,  we  have 
the  generator 


=  x  ~  -(2^x  +  uj2 x  )  ~  +  —L_ 

2  3x^  2  1  jx^  2  3x2 


(3.28) 


from  which  we  can  obtain  all  moments,  as  in  Section  II. 

Again  for  the  physical  noise  coefficient  system  (3.23),  we  cannot  obtain 
the  moments  directly.  However,  the  associated  Ito  system  (3.26)  possesses  the 


generator 


2  2-2 
(Y>_  „  3  r  0r  a\  .  2  ,  3  .  g*  r 

ic-  x2  -f(2cu-  2  )x  +  j  x  ]  ,x  +  2  a 

1  2  x. 


(3.29) 


from  which  all  moments  may  be  obtained. 

Finally,  it  is  immediately  seen  from  (3.19)  that  the  backward  operator  for 
the  system  (3.19),  which  contains  the  correction  terms  is 


j?-  T 


n  2  n  m  3g. ^ (x,t) 

:  [£  (x.t)+-j  I  I  g  (x,t)  -b- - 1  £ 

1=1  k=l  1=1  k  l 


l  n  I’m  I  2 

2,’.  ,  I  ,  8ir^X’t^(’js<X,t)°rs  TT*. 
i,j=l  J_r,s=l  J  I  i  j 


(3.30) 


For  a  comprehensive  study  of  the  sample  behavior  of  such  physical  noise  co¬ 
efficient  linear  systems  and  their  associated  Ito  systems  see  [3.8]. 
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In  closing,  we  mention  that  physical  noise  coefficient  systems  with  small 
parameters  can  also  be  related  to  associated  Ito  systems  and  their  Markov  process 
solutions.  These  techniques  are  based  upon  so-called  averaging  methods,  which 
will  be  discussed  in  detail  in  Section  V.  We  shall  now  turn  to  a  more  detailed 
discussion  of  the  physical  noise  coefficient  case. 
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IV.  Physical  Noise  Coefficient  Systems 

For  the  case  of  linear  systems  that  cannot  be  in  some  way  approximated  by 
differential  equations  with  white  noise  coefficients  or  approximated  by  other 
methods,  there  are  very  few  specific  results  available. 

This,  of  course,  is  due  to  the  fact  that  differential  equations  with  time 
varying  coefficients  cannot,  in  general,  be  solved  exactly.  Indeed,  unless  we 
can  solve  the  equations  exactly,  we  cannot  expect  to  obtain  moments,  nor 
probability  densities.  This  is  the  opposite  of  the  white  noise  coefficient  case, 
where  moments  can  be  obtained  exactly  even  though  the  differential  equations 
cannot  be  solved. 

The  question  that  we  must  first  consider  is,  what  types  of  ordinary  dif¬ 
ferential  equations  with  time  varying  coefficients  can  we  solve  exactly. 

The  general  equation  of  interest  is 


x ( t )  =  A(t)  x ( t ) , 


(4.1) 


where  x  is  an  n-vector,  and  A(t)  is  an  nxn  matrix  which  contains  elements 
that  are  randomly  time  varying. 

For  the  first  order  case 


x(t)  =  a(t)  x ( t ) , 


(4.2) 


we  can  write  the  solution  as 


a(s)ds 


x(t)  =  x  e 
o 


(4.3) 


whose  moment  properties  may  be  obtained  under  certain  conditions. 

This  depends  upon  the  assumptions  we  place  upon  the  coefficient  process 
(a(t) ,  t  €[0,°°) } . 

It  is  a  rather  important  point  that  in  general,  given  the  statistical 

properties  of  the  a-process,  (e.g.  the  joint  probability  densities)  the 

ft 

probability  densities  of  the  integrated  process  a(s)ds  cannot  be  obtained. 


-46- 


55 

•ft 


I 


'•s 


* 


§ 


,  * 

{  v 


.V 


% 

»  p 


! 


Thus,  even  for  the  first  order  physical  noise  coefficient  differential  equation, 
we  may  not  be  able  to  obtain  even  its  moments,  exactly. 


The  reason  here  is  that  we  must  evaluate 

ft 

a(s)ds 

i 

_  n‘ n  n 


a(s)ds 


E{x  (t) }  =  E{x  e 
o 


}  =  xnE{en'° 


}. 


(4.4) 


Without  knowledge  of  the  probability  density  for  the  integrated  a-process,  we 
cannot  explicitly  evaluate  these  expectations. 


Fortunately,  there  is  a  class  of  coefficient  processes  for  which  we  can 
determine  the  joint  probabilities  for  the  integrated  process.  This  class  is, 
of  course,  the  class  of  Gaussian  processes.  For  the  Gaussian  processes  we 
know  that  any  linear  operation  on  the  process  will  again  yield  a  Gaussian  process. 
This  is  the  only  general  class  for  which  we  can  make  such  a  statement.  Other 
processes  simply  do  not  allow  us  to  make  such  a  complete  statement.  Hence,  this 
is  another  case  in  which  the  Gaussian  assumption  allows  us  to  make  a  complete 
statement  about  the  solution  process. 


If  the  a-process  is  zero  mean  Gaussian  with  covariance  y(t  ,t  ),  it 

ft  a  l  z 


immediately  follows  that  the  integral, 


a(s)ds  is  also  a  zero  mean  Gaussian 


process  with  covariance  given  as 

t,  t 


ri. 


E  { j  a(s)ds 
•  o 


r  2  fl 

a(s)ds)  =  E  { 
o 


ds1  j  ds2  a(s1)a(s2)} 

o  *  o 


t, 

\  ' 


I  ds^  ds2  E{a(s^)a(s2) } 
^  o  '  O  L 


r"1 

Jo"1  -o 


ds,  j  ds?  ya(s1,s2)  . 


(4.5) 


Thus,  it  follows  that  the  integrated  a-process  is  zero  mean  Gaussian  with 
variance 


*  4 

o 

rt 

i 

dsi 

O 

ft 


ds2  Ya(srs2) 


(4.6) 
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and  the  desired  expectations  (4.4)  can  be  obtained  explicitly  as 

ft 

~  a(s)ds 

E{xZ(t)}  =  xn  E(en  lo 


} 


o 


,ny 


n 

=  x  e 
o 


—  CO  /2TTCT(t) 

T  a2(t) 


_ _ 1 _ 

2a2  (t)  , 

e  dy 


(4.7) 


where  a2  is  given  by  (4.6). 


For  the  non-Gaussian  case  one  might  consider  an  expansion  of  the 
exponential,  leading  us  to  the  result 

t 


E{e 


n  1  o 


a(s)ds 


nk  E{( 


k=0 


a(s)ds)  } 


k! 


(4.8) 


Since  the  moments  of  the  integrated  a-process  can  be  written  as 


E{( 


^  k. 

rc 

a(s)ds)  }  = 
o  J 

ds,  . . . 
o  1  J 

dsk  E(a(s1)...a(s  )}, 


(4.9) 


then  if  the  joint  moments  of  the  a-process  are  known,  we  can  in  principle 
write  (4.8)  in  a  series  form.  However,  summing  the  series  would  be  a  problem 
of  higher  order  of  magnitude  difficulty. 


Thus,  we  see  that  only  for  the  Gaussian  assumption  can  we  obtain  explicit 
results  for  even  the  simpHest  first  order  physical  noise  coefficient  system. 


What  can  be  said  for  higher  order  systems?  The  underlying  difficulty 

th 


here  is  simply  that  the  general  solution  to  the  n  order  time  varying  system 
(4.1)  cannot  be  written  in  a  closed  form.  One  is  tempted  to  write  the  solution 
as  the  matrix  exponential 

ft 

A(s)ds 


x(t)  =  eJo 


(4.10) 
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However,  as  is  well  known  (4.10)  can  be  the  solution  of  (4.1)  only  if 

ft 

a  condition  such  as  the  commutativity  of  the  matrices  A(t)  and  A(s)ds,  holds. 


Otherwise,  we  cannot  represent  the  solution  as  (4.10).  Of  course  this 
does  allow  us  to  study  systems  satisfying  the  equality 


A(t)  =  a(t)  A  , 
o 


(4.11) 


where  a(t)  is  a  scalar  random  process,  and  Aq  is  a  constant  matrix.  In  this 
case  A(t)  and  A(s)ds  =/ A  a(s)dsldo  commute,  so  that  the  solution  of  (4.1) 

VoJ“  ) 


with  (4.11)  is  given  by 


A  a(s)ds 

°  J 

x(t)  =  e  X 


(4.12) 


This  is  of  course  no  different  than  in  solving  for  the  solution  of  the 
constant  coefficient  matrix  case.  The  main  distinction  being  that  the 
characteristic  values  of  the  A  matrix  will  pick  up  factors  of  [  a(s)ds.  There¬ 


fore,  if  Aq,  for  example  can  be  diagonalized  by  the  matrix  X,  then  we  know 


X  1A  X  =  A 
o 


(4.13) 


where  A  is  the  diagonal  matrix  of  the  characteristic  values  of  Aq. 


Hence  (4.13)  yields 


X  A  a(s)ds  X  =  A  a(s)ds 
0  '  o  o 


a(s)ds  (^~^) 


O  \  f 


a(s)ds 


(4.14) 
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which  gives  us  the  explicit  solution 

.  re. 


x(t)  =  X  ^  e 


Jo 


a(s)ds 


X  x 


5 


r, 


(4.15) 


A 


i 

a 


I 


* 


a. 


a 


w> 
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It  follows  that  if  the  a-process  is  Gaussian,  then  as  in  the  first  order 
example,  we  can  obtain  all  of  the  desired  statistical  properties  of  the  solution 
vector  process  x(t). 

Can  we  say  anything  specific  for  more  general  noise  coefficient  linear 
systems? 

There  is  a  class  of  linear  random  coefficient  systems  for  which  we  can 
find  explicit  solutions.  In  order  to  describe  this  class,  we  first  appeal  to 
our  usual  first  order  case. 

Let  us  consider  the  first  order  equation  that  generalizes  (4.2),  which 
we  may  write  as. 


-  =  cx(t)  +  (l  d  a  (t))x(t)  , 
i=l 

where  c,  {d^-  are  given  constants  and  (a^(t)l  are  m-coeff icient  processes 
which  may  or  may  not  be  correlated,  may  or  may  not  be  Gaussian,  etc. 

We  can  write  the  solution,  with  x(0)  =  x  ,  as 

o 


(4.16) 


c  t  + 


x(t)  =  x  e 
o 


m  rt 

J,  diJ  ai(5) 

1=1  ;  o 


ds 


9 


.  -r. 


•  •  ^ 


(4.17) 


Having  this  explicit  solution,  we  can  therefore  obtain  statistical  properties  as 

ft 

desired  depending  upon  whether  the  joint  statistics  of  the  integrals  a.(s)ds 
can  be  obtained  (such  as  in  the  Gaussian  case). 


It  immediately  follows  that  if  the  matrices  C,  D_,,...,D  are  nxn  diagonal 

1  m 


matrices , 


then,  the  vector  equation 


=  C  x(t)  +  [l  a.(t)D.]  x(t) , 


may  be  written  as  the  n  first-order  equations 


dx  (t) 


=  c  x  (t)  +  [/,  d  a.(t)]x.(t),  j  =  1 ,2 , . .  .  ,n . 
J  J  i=l  1J  1  J 


The  solutions  to  (4.19)  may  be  written  as  in  the  form  of  (4.17). 

The  next  level  of  complexity  will  occur  when  the  C,  { D^}  matrices  are  all 
upper  triangular.  In  this  case,  the  solutions  may  be  obtained  sequentially. 


Thus  for 
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we  would  have 


m 


dx  (t) 

57 -  =  C„x^c>  +  tl  d^V01*^0 


dt 


n  n 


i=l 


dx  .  (t) 
n-1 

dt 


m 


'  VlVl(t|  +  c;.lxn(t>  +  t|_1di(n-l)a1(c>J  Vl(t)  + 


m 


[?  1din-l3i(t)]  xn(t) 
i=l 


id 


f 

I 

I 

r» 

k 


•jt 


etc. 


(4.21) 


We  see  that  there  is  a  sequence  of  first  order  equations,  all  solvable. 

The  first  equation  is  first  order  homogeneous,  all  others  are  first  order  non- 
horaogeneous.  Again,  if  all  the  [a^(t)]  are  Gaussian  processes,  then  the  statistical 
properties  of  the  solution  vector  process  x  are  completely  determined. 

The  next  step  in  the  generalization  is  to  consider  the  case  that  the 
matrices  C,{D^}  can  all  be  transformed  into  upper  triangular  matrices.  This  would 
allow  us  to  reduce  the  problem  to  the  form  we  have  just  considered.  This  problem 
can  be  expressed  in  a  concise  form  by  means  of  Lie  algebra  theory,  [4.1],  [4.2], 
[4.3],  [4.4]. 

We  briefly  describe  the  fundamental  ideas.  A  subspace  L  of  nxn  matrices 

is  called  a  Lie  algebra  if  for  all  C,  D  in  L,  the  commutator  product  [C,D]  =  CD-DC 

also  belongs  to  L.  Let  L(C,D,,...,D  )  denote  the  smallest  Lie  algebra  containing 

1  m 

the  matrices  C,  D^,...Dm.  This  is  usually  referred  to  as  the  Lie  algebra  generated 

bv  C,  D,,...D  .  One  defines  the  derived  series  of  the  algebra  as  follows, 
i  m 


i 


1 

I 


iti 


£ 


hi 


.*>**;«.*>, 
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L°  =  L 


L  =  [L,  L] 


Ln+1  =  [Ln, 


Ln] 


where  [L^,  ]  denotes  the  set  of  all  matrices  formed  by  [C^,D^]  =  C^D^-D^C^, 

,  „k  k  _  k 
where  C  ,  D  G  L  . 

The  Lie  algebra  is  said  to  be  solvable  if  there  is  an  integer  N  for 
N 

which  L  =  {0}.  Clearly,  an  Abelian  Lie  algebra  (where  all  matrices  are  pairwise 
commutative)  is  a  special  case  of  a  solvable  Lie  algebra,  since  C,  DSL  satisfies 
[C,  D]  =  {0}.  Therefore,  [L,  L]  =  {0}. 

How  do  these  results  impact  upon  our  problem  for  random  coefficient 
systems?  The  connection  with  our  discussions  above  is  contained  in  the 
following  lemma.  [4.5]. 

Lemma  -  A  matrix  Lie  algebra,  L,  is  solvable  if, and  only  if,  there  exists  a 
non-singular  matrix  P  such  that  PMP  ^  is  upper  triangular  for  any  M  £L- 

Thus,  we  see  that  if  the  nxn  matrices  C,  D.,...,D  generate- a  solvable 

1  m 

Lie  algebra,  L,  then  the  linear  system 


dx(t) 

dt 


in 

Cx(t)+  [£  a .  (t  )D  ]  x  (t ) , 
i=l  1 


(4.22) 


can  be  transformed  by  setting  y =  Px  where  P  is  guaranteed  to  exist  by  the  lemma 
to  yield 


dy(t)  _ 
dt 


PCP-1  y (t )  +  [l  a  (t) 
i=l 


m 


PD.P"1]  y (t)  , 


(4.23) 


in  which  each  matrix  PCP  ,  PD  P  ,  ...,PD  P  is  upper  triangular.  It  follows 

m 

that  the  solution  vector  y(t)  can  be  obtained  explicitly  by  sequential 
calculations.  Finally  the  x-process  is  obtained  by  x(t) =  P  ^y(t). 

Example  [4.2] 

Consider  the  second  order  system 


“Jt^  =  Dia1(t)  +  D2a2(t)]x(t), 


where 


C  =  c» 


2  -1 


1  0 


•D1  = 


1  0 


1  0 


.  °2  = 


1  -1 


0  0 


and  a^(t),  a2(t)  are  Gaussian  processes  possibly  correlated. 


It  is  possible  to_show  that  the  Lie  algebra  generated  by  (C,  D  ,  D  ) 
0  1  1  12 
is  solvable,  with  P  =  . 

1  -1 


We  immediately  find 


PCP  1  =  a 


1  1 


0  1 


,  PDXP 


1  1 


0  0 


,  pd2p 


0  0 


0  1 


which  gives  us  the  transformed  equation 


dv(t) 

dt 


‘  +  a1  (t) 


i+  a1(t) 


u+  a?(t) 
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Clearly,  (4.24)  can  be  explicitly  integrated,  giving  us  the  desired  statistics 
for  the  vector  solution  process  x=p  ^y. 

There  is  one  point  that  must  be  made  here  concerning  the  Lie  algebra 
approach.  Generally  speaking,  many  of  the  linear  equations  that  we  study  in 
mechanical  vibrations  do  not  generate  solvable  Lie  algebras. 

Thus,  the  classical  oscillator  equation  yields 


0  1 


0  0 


,  D  « 

-to2  -2?to  ;  1  0 


for  a  randomly  varying  stiffness  coefficient. 


x(t)  +  2^tox(t)  +  [to  +  a(t)]x(t)  =  0 


Even  in  this  simplest  case  therefore  we  cannot  use  the  Lie  algebra  approach. 
For  general  linear  systems,  one  would  have  to  test  whether  they  generate  solvable 
algebras.  This  is  not  easy  to  accomplish,  especially  for  structures  with  many 
components.  For  further  details  one  may  see  [4.5],  [4.6]. 

We  shall  finally  consider  what  can  be  stated  concerning  bounds  on  the 
second  moment  statistics  from  a  Lyapunov  function  point  of  view.  We  shall  follow 
the  derivation  of  Infante  [4.7]  [see  [4.8],  also].  For  deterministic  systems 
this  basic  idea  was  first  studied  by  Wintner  [4.9]. 

t  h 

Towards  this  end,  for  the  arbitrary  n  order  linear  system 


=  (A+ F(t))x(t)  , 


(4.25) 


We  define  the  quadratic  form 


V(x)  =  x  Px, 


(4.26) 


3 
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where  P  is  a  positive  definite  symmetric  matrix.  Clearly  if  P  is  the  identity, 
then 


V  (x)  =  ||  x  ||  ^  ,  ||  x  ||  =  (£x.V/2.  (4.27) 

We  evaluate  the  derivative  of  (4.26)  along  the  trajectories  of  (4.25), 
thus  obtaining. 


j-  V  (x(t) )  =  xT(t)  [(A+F(t))TP+P(A+F(t))]x(t) 


(4.28) 


Upon  dividing  by  V(x(t)),  we  will  obtain 


V(x(t)) 


dV(x(t))  _  x*(t)[(A  + 
dt 


F(t)  ) 
T 


P  +  P(A  +  F(t))]x(t) 
x  (t)  Px  (t) 


(4.29) 


Clearly  the  quotient  on  the  right  hand  side  of  (4.29)  is  quite  complicated. 
However,  by  the  properties  of  pencils  of  quadratic  forms  [4.10],  as  well  as  the 
min-max  properties  of  eigenvalues,  it  is  known  that 


X  .  (DP 
min 


< 


T 

x  Dx 
T 

x  Px 


< 


X  (DP 
max 


), 


(4.30) 


where  X  ,  X  are,  respectively,  the  minimum  and  maximum  eigenvalues 
in  in  max  ^  ^ 

of  the  matrix  DP  .  Note  since  P  is  positive  definite,  then  P  exists. 
Furthermore,  the  operator  associated  with  the  matrix  DP  ^  can  be  shown  to  be 
symmetric  [4.9],  hence  the  eigenvalues  of  DP  ^  must  be  real. 


that 


Upon  applying  the  inequalities  (4.30)  to  the  equality  (4.29),  we  find 


X  . 
mm 


(t) 


< 


1 

V(x(t ) ) 


dV(x(t) ) 
dt 


<  '  (t), 

—  max 


(4.31) 
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where 


Xn,-in  mav(C)  =  Kin  mflv  [  (A  +  F  (t)  )  T  +  P  (A  +  F  ( t  )  )  P_L  ]  . 
min  $  H13X  min  y  msx 

Upon  integrating  (4.31)  and  taking  exponentials  of  each  term,  we  have 


ft 


X  (s)ds 
min 


X  (s)ds 
max 


<  V(x(t ) )  <  e 


which  leads  to  the  expected  values 


E{< 


rt 


X  .  (s)ds 
mxn 


X  (s)ds 
max 


}  <  E{V(x(t) ) }  <  E{e 


}  . 


(4.32) 


For  the  left  inequality  of  (4.32),  we  may  apply  Jensen's  theorem  [4.11] 
which  states  that  if  g(x)  is  convex  (i.e.,  g"(x)  _>  0),  then  g(E{x})  <_  E{g(x)}. 

Thus  we  have  from  (4.32),  the  lower  bound 


ft 

E{  I  X  (s)ds } 

Jo  mln 

e  <  E(e 


fc 

^  o 


J  "min'8*'19 

J  r\ 


}  <  E{ V(x(t ) )  } 


(4.33) 


Therefore,  we  must  find 


E  i 


rc 

I 

'  o 


min 


(s)ds ) 


Ev  >  .  (s) }ds  . 

i  min 

‘  o 


For  the  oscillator 


x  +  2v,x  +  (a,  +f(t))x=0. 


we  have 


(4.34) 
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In  the  case  that  P=  I,  then  we  have  A  .  are  obtained  from,  the  eigen- 

min.max 

values  of  the  matrix, 


1-(uj2  +  f  (t) 


l-(u2  +  f(t) 


These  are  respectively 


co j  ±  i  c2<j2  +  (TM  +  W  . 


Hence,  the  lower  bound  on  the  E{  !(x  ||2  },  is  obtained  from  (4.33)  as 


Cwt—  E{/c2u2+  (f  (s)  +  to2-l )  ^  }  ds 


<  E{  II  x  II2  } 


By  the  Schwartz  inequality,  we  can  write 


E {  1 1 x il  }  <  E{  ixir> 


Thus,  it  follows  that 


(4.35) 


(4.36) 


E{/Ozu2  +  (f(s)  +  w2-l)2  )  -  E(c2'-2  +  (f(s) +u2-l)2}  7  ,  (4.37) 
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and  (4.36)  can  finally  be  written  as 


(i>-f;2u)2)t  -  E{  (f(s)  +u)2-l)2}ds 
e  <  E {  1 1  x  ( t )  1 1  Z}  . 


(4.38) 


I 


w 


£ 


.i 


Thus,  we  have  lower  bounds  that  can  be  evaluated  knowing  only  the  second 
moments  of  the  coefficient  process. 

Similar  results  can  be  obtained  for  upper  bounds  as  well  by  the  same 
concepts.  For  example,  see  [4.12]. 

This  section  contains  essentiallv  all  general  concents  that  allow 
explicit  solutions  for  the  moments  of  linear  systems  with  physical  noise 
coeff icients.  Any  other  exact  solutions  are  for  specific  cases.  In  general 
for  the  physical  noise  coefficient  case,  one  must  approximate.  This  is  the 
subject  of  the  next  section. 
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v.  Approximate  Methods 

Except  for  the  various  cases  and  examples  that  we  have  discussed  in  the 
previous  sections  for  which  exact  results  are  available,  all  other  linear  systems 
with  random  coefficients  must  be  looked  at  based  upon  some  approximate  method. 

These  methods  of  approximation  involve  expansions  in  small  parameters,  asymptotic 
convergences  in  terms  of  small  parameters  in  particular  the  various  stochastic 
averaging  techniques,  and  finally  the  so-called  closure  methods  in  which  expansions 
in  moments  or  cumulants  are  terminated  based  upon  some  assumed  property.  These  are 
the  main  topics  that  we  will  attempt  to  describe  in  this  section.  Since  exact 
results  are  few,  one  can  be  reasonably  comprehensive  in  a  survey.  But  approximate 
results  are  many,  therefore,  we  cannot  claim  to  be  comprehensive.  Instead,  we  try 
to  present  those  results  that  are  representative  of  the  specific  approach. 

V. 1  Expansion  Techniques 

As  we  stated  in  Sections  II,  IV,  when  the  coefficient  process  is  a  physical 
(non-white)  process,  we  cannot  separate  its  statistics  at  a  given  time  from  the 
response  statistics  at  a  given  time.  Therefore,  for  the  general  physical  noise 
we  cannot,  in  general,  even  determine  the  first  moment  exactly.  We  can  only 
approximate.  One  method  of  approximation  is  by  expansion.  The  literature 
contains  a  great  many  examples  of  this  approach. 

Basically  the  idea  is  as  follows. 

For  the  linear  system: 


=  (A+F(t))x(t),  (5.1) 

as  stated  before,  A  is  an  nxn  constant  matrix  and  F(t)  is  an  nxn  matrix 
whose  non-zero  elements  are  stochastic  processes.  We  can  write  (5.1)  as  the 
integral  equation,  with  initial  condition  x(t  ) =  x  , 


x  ( t ) 


Kt,t  ) x  + 
o  o 


ft 

t 

o 


$(t ,s)F(s)x(s)ds, 


(5.2) 


where  the  transition  matrix  0,  is  given  explicitly  as  the  deterministic  matrix 
exponential 


<t>(t,s) 


=  eA(t-s) 


(5 


where  (t-s)  is  a  scalar  parameter. 

A  typical  recursive  approximation  scheme  represents  the  solution  to  (5.2) 

as. 


$(t,t  )x 
o  o 


$(t,t  )x  +  |  v(t,s)F(s)x  (s)ds 
~  o 


c  )x  +  I 

o  o  Jc 


t>(t,t  )x  + 
o  o 


, s)F(s)x  - (s)ds 
n-X 


(5 


To  a  certain  extent  most  approximations  to  the  time  varying  linear  system  is  of 
this  form,  or  closely  related  to  it.  An  associated  procedure  is  to  consider  the 
expansion  in  e  of  the  system. 


dx(t ) 
dt 


(A  + eF(t))x(t)  , 


(5 


Upon  writing  the  associated  integral  equation  equivalent  to  (5.2), 
we  have 

t 

x(t)  =  ?(Cit  )x  +  e  I  $(t,s)F(s)x(s)ds.  (5 


-61- 


Now  expanding  Che  solution  x(t)  in  powers  of  e,  we  obtain  the  representation 


:(t)  =  l  e  x,  (t) 


(5.7) 


Substituting  (5.7)  into  (5.6)  and  identifying  the  coefficients  of  the  powers 
of  c ,  we  obtain  the  sequence  of  terms  x^(t),  as 


xQ(t)  =  Ht,to)xo 


x.(t)  =  *(t,s)F(s)x  (s)ds 

1  it  0 

o 


x2(t)  =  j  $(t,s)F(s)x^(s)ds 
.  Co 


x  (t)  =  <t(t  ,s)F(s)x  T  (s)ds 
J  t  n-1 

.  o 


(5.8) 


Both  of  these  ideas  are  quite  classical.  For  example,  one  can  find  the 
representation  (5.7)  and  expansion  (5.8)  in  the  works  of  Kryloff  and  Bogoliubov 
[5.1].  For  a  discussion  see  [5.2].  Indeed,  this  procedure  is  attributed  to  Poisson 
and  was  studied  further  by  Poincare  in  1892.  Naturally  convergence  is  the  problem 
of  significance,  especially  if  e  is  not  small.  For  e  =  1,  this  has  been  called 
Adomian's  method  in  the  stochastic  literature.  See  for  example  [5.3]  and  [5.4]. 
Given  such  a  classical  approach,  it  appears  somewhat  questionable  to  attribute 
the  procedure  to  anyone  currently  active. 

Of  course  the  objective  of  these  expansions  is  to  obtain  the  statistics 
of  the  solution  process.  Following  the  notation  in  [5.3]  we  can  write  (5.2)  as 


v  v  *v;\vy 


m 


x(t)  -  1  N(t , s)x(s) ds  =  G(t ,  t  ), 
J  t-  ° 


where 


N(t,s)  =  -i>(t ,  s)  F(s) 

G(t ,  t  )  =  $(t,t  ) x 
o  oo 


Defining  the  iterated  kernels  N^(t,s)  as 


N^t.s)  =  N ( t , s) 


Nk(t,s)  =  j  N(t,x)Nk_1(r,s)dT 
Co 


k  =  2,3,.  .. 


we  define  the  resolvent  kernel  of  N(t,s)  by  the  Neuman  series 


T(t,s)  =  l  N  (t,s) 


This  allows  one  to  write  the  solution  to  (5.2)  as 


x ( t )  =  4> ( t , t  ) x  +  |  r(t,s)<i>(s,t  )x  ds. 

o  o  J  o  o 

o 


In  this  case,  the  mean  vector  E{x(t)},  may  be  obtained  as 


E{x(t)}  =  ‘,(t,t  )E{x  }  +  E{  F(t  ,s)  (s,t  )E{x  ids, 
o  o  o  o 


assuming  that  the  statistics  of  F(t),  and  of  xq  are  indepei dent  of  one  another, 

and,  furthermore,  that  the  A  matrix  is  deterministic.  If  the  A  matrix  contains 

random  constants,  then  the  terms  E{i(t,t  )},  E{$(t,s)}  would  appear  in  (5.13). 

o 

Although  these  procedures  are  straightforward  applications  of  established 
solution  expansion  techniques,  obtaining  the  desired  statistics  is  not  made 
truly  simpler  due  to  the  complexities  of  the  infinite  number  of  terms  to  be 
considered.  Expansions  of  moments  using  such  ideas  have  been  studied  by  many 
investigators.  In  the  modern  development  of  the  topic,  one  can  go  back  to  the 
works  of  Samuels  and  Eringen  [5.5],  [5.6]  who  studied  stability  of  the  second 
moments  of  systems  of  the  form  (5.1)  by  such  expansion  techniques,  (see  also 
[5.3],  [5.7]. 

For  a  general  development  of  this  topic,  one  might  look  at  [5.8],  Notice 
that  there  is  nothing  fundamentally  probabilistic  about  the  expansion  methods 
above.  The  so-called  Neuman  series  for  integral  equations,  as  well  as  the 
perturbation  expansions  were  developed  for  classical  deterministic  equations. 

There  are  a  number  of  papers  which  deal  directly  with  approximations  to 
the  joint  density  functions  via  a  sequence  of  integro-dif ferential  equations, 
that  involve  certain  smoothing  operations. 

The  first  paper  along  these  lines  apparently  is  due  to  Kryloff  and  Bogoliubov 
in  1939  [5.9].  Subsequent  studies  appeared  in  [5.10]  and  in  [5.11].  Also  consult 
[5.12]  for  further  referenced  details. 

V.2  Hierarchy  Equations  and  Closure 

Under  certain  assumptions,  it  is  possible  to  obtain  a  set  of  equations 
that  will  yield  the  moments  of  the  solutions  of  systems  with  randomly  varying 
parameters.  However,  these  equations  may  depend  upon  an  infinite  number  of 
variables  that  render  them  impossible  to  solve  by  finite  methods.  This  is 
typical  of  the  type  of  problem  one  meets  with  non-linear  systems  in  which 
lower  order  approximants  may  be  a  function  of  higher  order  approximants  so  that 
the  equations  are  not  closed.  In  that  case,  assumptions  are  put  on  the  system 
(which  may  or  may  not  be  justified)  which  will  reduce  the  infinite  heirarchv  of 
equations  to  a  finite  set  that  can  be  solved.  After  solution  of  the  finite  set, 
naturally  the  errors  in  the  approximation  must  be  studied.  At  this  time,  however, 
there  is  no  general  statement  that  can  be  made  regarding  the  errors  introduced. 


We  illustrate  this  problem  for  the  case  of  the  oscillator  with  a  randomly  varying 
spring  coefficient.  We  assume  that  the  random  coefficient  is  generated  by  filtering 
a  white  noise.  Thus,  it  is  a  Markov  process,  or  more  generally,  it  is  the  component 
of  a  Markov  process.  One  of  the  earliest  papers  to  treat  the  problem  in  this  way 
was  [5.13].  Thus,  consider  the  oscillator. 


x(t)  +  2cux(t)  +  [>j  +  y(t)]x(t)  =  0, 


(5.14) 


where  v(t)  is  generated  by  the  Ito  equation 


dy ( t )  =~by(t )dt +  dB(t) , 


where  the  B-process  is  the  usual  Brownian  motion.  We  recognize  that  the 
y-process  is  the  Ornstein-U'nlenbeck  process. 

Now  since  y  is  generated  by  an  Ito  equation,  then  it  follows  that  the 
vector  (x,  x,  y)  is  a  Markov  process  generated  by  the  Ito  system 


dx^(t)  =  x^ (t)dt 


dx£  ( t )  =-[2Cojx9  (t)  +  +  y(t))x1(t)  ]dt 


dt(t)  =~Sy(t)dt  +  dB(t), 


(5.15) 


(5.16) 


Following  along  the  concepts  developed  in  Section  II,  we  can  write  the  generator 
of  the  system  (5.16)  as 


$,-*2  sf  -[2wx2+(„2  +  y)x1] 


3  ,  J-  J- 

y  3y  2  3y,  ’ 


(5.17) 


and  obtain  the  moment  equations  via  the  Dynkin  formula,  (2.22),  as  in 
Section  II. 


I 


Upon  denoting  the  n  moment  as 


E{x^(t)x2(t)y  (t)}  =  nij ^  ( t )  ,  where  j+k+  1=  m 


we  shall  find  the  first  moment  equations  to  be 


mioo  moio 


moio  ”  mioo  2c'jmoio~mioi 


m001  Bm001  ‘ 


Here  we  already  see  that  the  m  ^  equation  contains  a  second  moment. 


mioi  =  E'x1(t)y(t)> 


(5.ia; 


(5.19) 


For  the  general  case,  we  shall  have 


rajkc  =  j  mj -1 , k  +  1 , £  "2Cukrajk£  "W  +  1 , k-1 , i  +  l,k-l,i  +  1  (5.20) 


-3*mjk£  +  -  l(l- l)m.k^_2 


In  (5.20),  we  see  the  term  ra  ,  which  is  j  +  k+?  +  l  moment,  one  order 

J  «  1  ^  K  J.  *  '■>  '  1 

higher  than  m  . 

J  iCa' 

Hence,  the  set  of  equations  m  for  j  +  k+£  <  n  is  not  a  closed  set.  It  is, 

j  kx 

instead  an  infinite  hierarchy  which  we  can  express  generically  as 


m  G  (m,,!!!-,...,!!!  ,m  ,  ,  )  , 
n  niz  nn+i 


(5.21) 


always  depending  upon  one  order  higher. 


rtf 


P  1  'V  I 


s 

I 

I 
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This  n,  of  course,  is  not  a  surprise  since  the  system  (5.16)  is  non-linear 
in  the  variables  (x^jX^.y)  through  the  x^y  term  in  the  second  equation. 

The  objective  of  a  closure  procedure,  would  be  to  invoke  a  property  or 
assumption  that  would  allow  us  to  write  the  approximate  equations 


B 

E 

B 

5 

I 


m  =  G  (m  ,m  , . . . ,ra  ).  (5.22) 

n  n  1  2  n 

There  is  no  general  theory  that  allows  one  to  measure  the  errors  made  with  such 
closure  procedures.  However,  there  have  been  many  attempts  to  study  this 
problem.  In  the  study  of  turbulence,  [5.14]  contains  investigations  of  closure. 
Closure  ideas  were  applied  to  stochastic  eigenvalue  problems  in  [5.15].  Lie 
algebra  ideas  were  applied  in  [5.16]  to  obtain  closed  first  and  second  moment 
equations . 

An  investigation  of  error  in  the  approximation,  leads  to  results  for  one¬ 
dimensional  systems  in  [5.17]  and  [5.18]. 

Two  classes  of  closure  concepts  naturally  arise.  They  are  Gaussian 
closure  methods  and  non-Gaussian  closure  methods. 

The  Gaussian  closure  idea  is  quite  clear.  Since  all  moments  of  a 
Gaussian  process  may  be  expressed  in  terms  of  the  first  and  second  moments, 
then  the  hierarchy  is  essentially  closed  by  writing  the  equations  for  the 
first  two  moments.  Any  higher  order  moments  that  appear  in  the  first  two 
equations  are  then  reduced  to  the  appropriate  expressions  in  the  first  and 
second  moments.  Thus,  for  example,  for  the  equation  (5.21)  Gaussian  closure 
reduces  it  to, 

m  =  H  (m- ,  m,)  (5.23) 

n  n  l  / 

for  some  new  function  H  ,  and  therefore 


m^  (m^ ,  m^) 


L  =  H  (m^,  m?) 


becomes  a  closed  set  that  can  be  solved  in  all  cases  for  linear  systems. 


W- fcit  NtaV . 


L  Vv  ■-.* 


•  kV  »  f  f  7  »  )  i 


(5.24) 
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Of  course  the  problem  that  exists  here  is  that  systems  with  randomly  varying 
coefficients  do  not  in  general  lead  to  Gaussian  solutions.  The  measure  of 
error  from  the  Gaussian  assumption  may  be  large.  This  in  turn  may  generate 
large  errors  from  the  time  moments.  The  reference  [5.19]  contains  further 
discussion  and  references  of  these  ideas  as  they  relate  to  structural  applications. 

Non-Gaussian  closure  techniques  assume  that  all  moments  of  order  r  >  n  can 
be  expressed  in  terms  of  lower  order  moments  r  £  n.  One  may  also  express  the 
non-Gaussian  closure  assumptions  in  terms  of  cumulants.  See  [5.19]  for  a 
discussion,  and  [5.20]  for  a  recent  application.  Cumulant  approximations 
have  been  used  in  studies  of  turbulence  [5.21],  and  many  other  stochastic 
continuing  problems  [5.22].  The  basic  idea  comes  from  the  characteristic 
function  of  the  underlying  probability  density  for  the  solution  process. 

To  that  end,  let  us  recall  that  the  characteristic  function  of  the 
process  x(t),  is  given  as 


(u,t)  -  EU1"*^}, 


(5.2! 


which  we  can  write  in  a  Taylor  expansion  about  the  origin  as 


r(u,t)  =  l  mk(t), 


(5.26 


k=0 


where  mk(t)  =  E{x  (t)},  assuming  all  moments  exist.  (Otherwise  (5.26) 


becomes  a  finite  series  with  a  remainder  term). 


An  equally  important  expansion  is  that  of  log  ?^(u,t),  which  we  shall 


write  as 


log 


:  (u,t) 

X 


GO 

\ 

L 

k=0 


(i  u  ) 
k! 


\(t) 


(5.27 


The  (t); 
x-process . 


are  referred  to  as  the  cumulants,  or  semi-invariants  of  the 

The  connection  between  m,  and  >  are  obtained  as  follows  [5.23], 

k  k 


K 


2  2 
1  =  ° 


A ^  -  3  m^m^  +  2m^ 


2  2  4 

A.  =  m  -  3  m.  -  4  m.m.  +  12  m.m„  -6m 
44  2  13  12  1 


(5.28) 


Inversely,  we  can  also  write 


A2  +  \ 


m3  *  A3  +  ^1A2  +  A1 


m.  =  A.  +  3A_  +  4A  A  +  6A^A  +  A^ 
4  4  2  13  12  1 


The  cumulants  are  also  a  measure  of  the  "distance"  away  from  Gaussian. 
For  example,  if  m^=0,  we  have  from  (5.29), 


(5.29) 


m2  =  A2<‘=  °  ^ 


“3  '  X3 


"4  '  31 2  +  '4 


"5  '  10X2-'3  +>5 


(5.30) 


But,  for  Gaussian  processes  one  finds 


=  x2^=  a  ^ 


=  0 


m  =  3A 
4  2. 


=  0 


The  way  in  which  the  so-called  cumulant  neglect  procedure  is  applied,  is  simply 
to  assume  that  for  some  K,  such  that  for  k  >  K,  we  set  A^=0. 

Since  A^  is  a  polynomial  in  m^,  m^  , . . .  ,m  ,  it  follows  that  this  will 
essentially  close  the  moment  equations.  Thus,  for  our  example  (5.14),  the 
second  moment  equations  are 


m200  2mllO 


m020  =  ^uino20  2w  m110  2mlll 


m110  m020  2?U)in110  m100  ~2m201 


m101  *  m011  -Smioi 


m011  2?wm011  “  "'lOl  ”  m102  m011 


m002  =  “2S  m022  +  a 


One  can  use  the  formula  for  the  characteristic  function  in  terms  of  the 
cumulants 


i(u  x  (t) +  u  x  (t) +  u  y(t)) 
<?(Ul,u7,u3,t)  =  E{e  J  } 


r  prs  ,p+r+s  p  r  s 

■  “O  l  .  1  U1  "2  "3 

p,r,s  =  0 


(5.33) 


I 


and  the  fact  that 


■1 

t  J 


.j+k+i  3j  3k  dl 

1  mjk«.  3U]j  3u2k  3u3i  ?(u1»u2*u3;t) 


to  obtain 


m100  X100 


m200  X200 


“ill  Xlll  +  xonxioo  +  xoioxioi+ xooixno+ xoioxioo 


Therefore,  if  one  assumes  that  the  third  order  cumulant  A  is  zero  the 
relation  (5.35)  would  yield 


(5.34) 


(5.35) 


mlll  =  ra011ml00  +  m010m!01  +  m001rall0  +  raoiomioo 


(5.36) 


which  relates  the  third  moment  m  ^  Co  first  and  second  order  moments. 
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In  the  same  fashion  we  could  find  similar  relations  for  m201  an<*  ml02’ 

Thus,  (5.32)  along  with  (5.19)  would  become  a  closed  set  of  non-linear  equations 
to  solve  for  the  moments.  This,  however,  would  require  numerical  methods  for  the 
solution.  Explicit  examples  of  these  concepts  can  be  found  in  [5.19].  For  further 
discussion  of  these  concepts  one  can  consult  [5.24]. 


V. 3  Methods  of  Averaging 

The  method  of  averaging  for  classical  deterministic  equations  has  its  origins 
in  the  works  of  Van  der  Pol  [5.25],  who  studied  the  non-linear  equation 


X  ( t )  +  Uj2x(t)  =  E  f  (x(t),x(t)). 


(5.37) 


where  e >  0  is  a  small  parameter. 

He  assumed  a  solution  in  the,  now  classical,  form 


x(t)  =  a(t)  cos(ut  +  f(t)). 


(5.38) 


where  a,  $  satisfy  the  equations 


a(t)  *  e  fj.  (a (t) ) 


v(t)  »  tf,  (a(t) ) 


(5.39) 


For  small  e,  they  are,  clearly,  slowly  varying  and  become  constants  as 


The  idea  was  applied  by  a  great  many  experts  in  non-linear  oscillations 
such  as  Mandelstam,  Papalexi,  Krylov  and  Bogoliubov,  Minorski  as  well  as  others. 
However,  the  method  of  averaging  was  not  placed  upon  a  firm  mathematical  foundation 
until  the  major  studies  of  Bogoliubov  [5.26]. 

There  are  many  interpretations  of  the  method  of  averaging  especially,  as  a 
result  of  the  concept  playing  an  increasingly  important  role  in  applications  to 
stochastic  problems.  For  a  general  survey  of  the  method  of  averaging  for 
deterministic  systems,  see  [5.27],  for  stochastic  systems  see  [5.19]. 


mmsm 


:  i*sW 


The  classic  result  of  Bogoliubov  is  as  follows:  Let  a  system  be  written 
in  "standard  form",  x, f  are  n-vectors 

x(t)  =  ef  (x(t),t),  x(0)  =  x  ,  (5 

o 

where  f  satisfies  boundedness  and  uniform  Lipschitz  conditions  in  x  for  some 
region  DCRn  (R  is  n-dimensional  euclidean  space).  Furthermore,  let  the 
limit  (time  average) 

i  rT 

f(y)  -  lim  —  f (y ,t)dt  (5 

Tt“  '  -*o 

exist  uniformly  for  y£D.  Then,  given  any  n  >  0,  and  an  arbitrarily  large 
T,  there  exists  an  £q,  such  that  for  0<e<£o>  the  solution  y  of  the  averaged 
equation, 

y(t)  =  e  f  (y) ,  y(0)  =  xq,  (5 

will  satisfy, 

j|  y(t)-x(t)  ||  ''  n  ,  (5 

for  t  (0,  T/s  ] . 

Of  even  more  interest  is  the  case  that  f(y,t)  is  almost  periodic  in  t, 
uniformly  for  y€D. 

In  that  case  it  is  a  classical  fact  that  the  limit  (5.41)  exists.  Further 
more,  Bogoliubov  established  that  the  inequality  (5.43)  holds  for  tf(0,"’). 
Therefore,  the  averaged  solution  y(t)  becomes  a  uniformly  close  approximation  to 
the  true  solution  x( t)  on  the  entire  time  domain  (0,°). 


This  is  an  especially  important  result,  since  it  allows  us  to  determine 
qualitative  properties  such  as  stability  and  boundedness  of  the  solutions  of  the 
original  system  (5.40),  in  terms  of  the  averaged  constant  coefficient  system 
(5.42). 

This  is  especially  straightforward  for  linear  systems  of  the  form 


x(t)  =  e  A  (t)x(t) .  (5 

In  order  to  apply  these  ideas,  a  linear  system  must  be  put  in  the  standard 
form  (5.44).  We  can  see  how  this  may  be  done  through  the  following  illustrative 
example. 

Example 

Consider  the  simple  oscillator  with  periodic  coefficient, 


^  -x-i— ^  +  (1  +  cosut)x(t)  =  0.  (5 

dt 

we  may  write  this  Mathieu  type  equation  as 

dx^(t) 
dt 

dx  (t) 

— - -  =  -(1  +  coswt)x(t).  (5. 

dt 

Let  us  assume  that  w  is  very  large  as  compared  to  unity.  Thus, 
coswt  is  a  very  rapidly  oscillating  coefficient  terra. 

In  particular  let  us  assume  that 


=  x2(t) 


0) 

•w  =  —  .  where  0  <  c . 


'  <  <  1 . 


(5 
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Therefore,  by  Bogliubov's  result  the  solution  vector  y(t)  defined  by 


(  dy  (x) 

Jr  “  y2(lJ 

■ 

dy?(x) 

^  — dr — ”“y1(T>»  y(°)  =  x(o)» 


(5.51) 


will  remain  "close"  to  the  solution  of  (5.46)  for  all  time,  t  £.  [0,°°). 

The  question  now  arises,  is  there  a  version  of  the  averaging  method  that 
is  possessed  by  systems 


x(t)  =  f(x(t),  n (t) ) ,  (5.52) 

where  n(t)  is  a  stochastic  process  with  suitable  properties. 

The  motivation  for  asking  this  question  is  quite  clear.  It  is  a 
fundamental  fact  that  if  the  elements  a„(t)  of  the  stochastic  matrix  A(t) 
are  stationary  stochastic  processes,  then  the  time  average 

1  fT 

lim  -  A(t)dt  (5.53) 

T^  •»  o 

exists  with  probability  one. 

Furthermore,  if  the  elements  a  (t)  are  ergodic  processes,  then  the 
limit  (5.53)  is  equal  to  E{A(t)}  with  probability  one.  As  a  result  of  this 
fact,  it  follows  that  there  should  be  an  averaging  result  for  stochastic 
systems. 


4, 

I 


V  v* 


We  can  illustrate  this  for  the  case  of  the  inverted  pendulum  [see  Fig.  l], 
with  base  excitation  n(t).  The  linearized  equation  is 

8  +  A  0-  \  (n(t)  +  g)  9  =  0  (5.54) 

m  i 

where  m,c,ii  are  mass,  damping  factor,  length,  and  g  is  the  gravitational 
constant . 

2L 

In  terms  of  canonical  variables,  8,  =  —  ,  where  L  =  T-V  is  the 

26 

Lagrangian  and  i  is  the  generalized  momentum,  we  obtain  the  equations 

?  =-m[mi)  ^ )  +  n(t)S  ]  [n(t)  +  2  i  —  ]  +  mgi  0 

m 

9  =  (mJ.“)  ■*"  ?  +  1  ^  n(t)9  (5.55) 

as  a  phase  space  form  of  (5.54). 

Now  let  us  assume  that  the  random  base  excitation  function  n(t)  is  of 
the  form, 

n(t)  =  cw (c  *t)  ,  (5.56) 

where  c  >  0  is  a  small  parameter. 

We  see  from  (5.56)  that  the  s  coefficient  implies  a  small  variance  for 
n(t)  if  the  variance  of  w  is  fixed.  Furthermore,  the  time  shift  t  -*•  e  ^t  in  w 
has  the  effect  of  shifting  the  average  power  to  higher  frequency  ranges. 

Upon  setting  i  \  ;  t  in  (5.56),  (5.55),  we  will  obtain,  with  'V 

denoting  — 
dr 
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-em[(m£)  +  w" (t)6]  [w' (t)  +  2  —  ]  +  emg£  9 

m 


ei  ^[(m£)  *")>  + w'(x)0] , 


(5.57) 


which  is  in  the  standard  form  to  apply  the  theorem  of  Bogoliubov. 
We  assume  that  w'(x)  is  stationary  ergodic  with 


E{ w ' (t ) }  =  0,  E{  (w ' (t  ) ) 2 }  =  a2. 


(5.58) 


Upon  writing  the  system  (5.57)  in  the  form  (5.44) ,  we  find  that 
the  A-matrix  is  simply 


A(t)  = 


-  i  [w'(t)  +  2  *£■] 
£  m 


-m  [w ' (t) 2  +  w'(t)  2  ~  -g.H] 


T  w"(t) 

X 


(5.59) 


whose  time  average,  by  the  assumed  ergodic  assumption  on  the  w  (x)  process,  is 
obtained  from  (5.58)  as, 


m  [a  -g i] 


(5.60) 
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Therefore,  Bogoliubov's  theorem  implies  that  the  sample  solutions  of  the 
original  system  (5.57)  will  be  close  to  the  solution  of  the  constant  coefficient 
deterministic  system  for  time  intervals  of  the  order  0(l/e).  We  note  an 
immediate  result  of  (5.60).  If  a2  >  gl,  then  A  is  a  stability  matrix.  Thus 
all  sample  solutions  of  (5.57)  should  decay  exponentially,  aL  least  for  time 
intervals  of  order  0(l/e). 

Simulation  studies  for  the  inverted  pendulum  with  stochastic  base  motion 
[5.28]  have  suggested  that  the  solutions  might  be  close  for  all  time  if  the 
solution  process  is  bounded  and  the  power  is  concentrated  in  the  higher  frequencies. 
A  proof  for  a  specific  form  of  base  excitation  that  is  bounded  was  given  by 
Bogdanoff  in  1962  [5.29].  He  considered  excitations  of  the  form 


n(t)  =  y  a.  cos(j.t+  ^.), 


(5.61) 


where  [a^i,  are  given  constants,  all  a^  are  "small"  in  magnitude, 

w.  are  "large"  and  all  [ -j  -cu.l  >K  for  some  constant  K. 
l  &  1  i  J 

Finally,  the  [v  }  are  independent  random  variables  all  uniformly 
distributed  on  [0,  211].  The  functional  form,  with  these  assumptions  does 
fit  within  the  scope  of  Bogoliubov's  theorem  for  almost  periodic  excitations. 
Bogdanoff's  proof  was  based  upon  1 inearizat ions  and  comparing  of  small  terms. 
Experimental  results  [5.30]  corroborated  the  theory  of  [5.29]  excellently,  even 
to  the  loss  of  stability  when  two  distinct  frequencies  become  close. 

A  generalization  of  the  results  in  [5.29]  was  presented  in  [5.31]  for 
linear  systems  in  the  standard  form  (5.44),  where  the  ergodic  matrix  A(t) 
satisfies. 


(a)  Sup  ||  A ( t )  ||  <  M, 


(b)  Sup  :  [A(t)-P]dt  '|  <  M0 , 

t  J  o 


(5.62) 


where  P  =  ErA(t) ; . 


It  was  established  that  if  P  is  a  stability  matrix,  then  the  solution  to 


y  =  £Py,  y(0)  =  x(0)  , 


satisfies 


lin  Sup  II  x(t)-y(t)  ||  =  0  , 

elO  t  e[0,«) 


with  probability  one,  where  x(t)  is  a  sample  solution  to  the  stochastic 
linear  system  (5.44). 

Unfortunately,  this  result  requires  the  rather  strong  assumption 
(5.62)(b).  Although,  this  is  known  to  hold  for  almost  periodic  matrices, 
A(t),  we  cannot  state  any  general  statistical  properties  that  will  guarantee 
this  assumption  to  hold  in  the  stochastic  case.  These  assumptions  were  also 
used  in  [5.32].  For  more  recent  results  see  [5.33]. 

We  should  mention  at  this  point  that  even  if  the  time  average  of  the 
stochastic  coefficient  matrix  A(t)  exists,  we  cannot  expect  the  solution 
process  to  be  close  to  the  solution  of  the  ensemble  averaged  system,  if 
there  are  no  further  assumptions  such  as  ergodicity.  This  can  be  seen  by  the 
following  simple  example. 

Example .  Consider  the  first  order  scalar  equation 


x  =  tbx,  x(0)  =  1 , 


where  b  is  a  bounded  constant  random  variable  with  E(b}  =  p<  0. 
Thus,  the  ensemble  averaged  system  is 


y  =  f;py,  y(0)  =  l . 


mm 


■ '  •■VSVr/Vv  ■; ,  1. 


1w*\v/r/fwnwv  ^ 


How  do  the  solutions,  x(t),  y(t)  compare?  In  particular,  what  can  be  said 
about 


lin  |x(t)  -  y(t)  |  =  lin  |  eebt  -  e  “  ^  t  \  . 
c4-0  e-f-0 

It  can  be  shown  that  for  a,3>0,  a^8,  we  have  by  simple  calculus. 


_ o_  a 

—Qt  ft  tt  — 8  „  o~8 

Sup  | e  -  e  |  =  ! (f)  -  (-)  |  >  0, 

Ct  UL 

c 


which  occurs  at  t  =  — —  (loga  -  log8)  >  0. 

For  our  case  (5.65),  the  e  cancels  out,  and  we  will  obtain  for  b <  0, 


b-P  b-P 

_  ,  sbt  ept,  |  /Pv  ■  ■  ,P,  . 

Sup  !  e  -  e  I  =  I  (t-)  -  (t-)  !  >  0 

t  b  b 


independent  of  e ! 

Hence,  the  solutions  cannot  become  close  as  t  1  0,  assuming  only  that 
time  averages  exist! 

For  the  case  of  the  white  noise  excited  system,  we  can  write  an 
averaging  result  that  relates  the  stochastic  response  characteristics  to  the 
associated  averaged  deterministic  system,  [5.34]. 

In  particular,  for  an  Ito  system  of  the  form. 


dx  =  c  (Axdt  +  F(x)dB ) ,  x(0)=x 

o 

where  A  is  a  constant  or  periodic  matrix,  F(x)  is  a  matrix  whose  elements 
are  linear  functions  of  the  vector  x,  and  B  denotes  a  vector  of  Browmian 
motions,  it  follows  that  the  averaged  system 
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y  =  eAy ,  y(0)  =  x 


(5.69) 


satisf ies 


lim  Prob  {  Sup  J|y(t)-x(t)||>_6}  =  0 
e+0  t  £  [  0 ,  °°) 


(5.70) 


for  any  5  >  0. 

In  this  case,  as  before,  A  denotes  the  time  averaged  matrix  defined 
by  (5.53). 

The  effective  importance  of  this  result,  is  that  if  the  structural 
excitations  are  wide  band  Gaussian  then  we  can  treat  the  system  (5.68)  for 
small  e,  as  the  deterministic  averaged  system  (5.69). 

Furthermore,  the  response  characteristics  of  the  true  system  (5.68) 
will  remain  close  to  the  solution  of  (5.69)  for  the  entire  time  interval 
[0,”) • 

The  question  of  controlling  the  true  system  through  the  averaged  model  also 
becomes  possible.  The  control  problem  is  discussed  in  [5.34]  as  well. 

Many  averaging  methods  that  have  been  developed  for  the  study  of  stochastic 
systems  have  not  taken  the  view  that  we  have  discussed  above.  The  basic 
philosophy  that  we  have  discussed  above  is  essentially  to  replace  the  true  system 
with  random  excitations  by  a  system  that  essentially  averages  out  the  random 
fluctuations  leaving,  in  a  sense,  the  mean  system  which  is  deterministic.  If 
the  true  and  the  mean  ^sterns  are  "close"  for  all  time,  then  we  can  restrict  our 
study  to  the  simpler  mean  system  response,  knowing  that  it  will  be  very  close  to 
the  true  system  dynamics.  On  the  other  hand,  the  philosophy  that  drives  the  work 
of  Khazrainski  [5-35],  [5-36],  [5-37],  Stratonovich  [5-38]  and  Papanicolaou  [5.39] 
is  to  replace  the  original  random  coefficient  system  by  one  that  is  simpler 
through  an  averaging  procedure.  Generally  speaking  their  averaged  systems  are 
close  to  the  true  systems  in  the  sense  that  the  probability  distributions 
generated  by  the  averaged  system  are  close  to  the  probability  distributions 
generated  by  the  true  system.  This  is  referred  to  in  probability  theory  as 
weak  convergence. 


It  differs  from  the  results  we  have  described  above,  in  which  we  are 
concerned  with  the  actual  response  of  the  system  being  close  to  the  response 
of  the  averaged  model.  This  is  called  strong  convergence,  or  almost  sure  sample 
convergence*  which  will  imply  that  the  probability  distributions  will  be  close. 
The  important  difference  is  that  from  the  structural  engineering  point  of  view 
we  are  concerned  primarily  with  the  actual  structural  response  of  the  systems. 
However,  the  weak  convergence  point  of  view  does  allow  one  to  obtain  probability 
distributions  for  averaged  systems,  when  it  may  not  be  possible  to  determine  the 
distributions  for  the  true  randomly  excited  system. 

To  a  great  extent  this  approach  was  motivated  by  [5-35],  where  an 
averaging  result  was  established  for  partial  differential  equations.  In 
particular,  we  consider 


3u  _  . 

—  =  eL(x,t)u, 


(5.71) 


where  L  is  an  elliptic  or  parabolic  second  order  differential  operat 
with  sufficient  regularity  conditions  relative  to  the  (x,t)  variables. 


L  (x)  =  lim  i  J  L(x , t ) dt , 
°  T+»  T  J0 


(5.72) 


be  the  time  averaged  operator. 

Then,  if  v(x,t)  is  the  solution  to 


3v 

eL0(x)v’ 


(5.73) 
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for  the  boundary  condition 


u(x,T/e)  =  v  (x,T/e)  =  f(x) 


(5.74) 


it  follows  that. 


lira  Sup  |u(x,t  -v  (x,t)  I  =  0 

e  +  0  ( x , t )  €.1^  [0,T/e] 


(5.75) 


The  importance  of  this  theorem  for  stochastically  perturbed  systems  is  the 
following;  for  the  Ito  system 


dx = m(x,t)dt + o(x,t)dB,  a  =  (o„  ) 


(5.76) 


where  x,m  are  n-vectors,  a  is  an  nxn  matrix  and  B  is  an  n-vector  of 


Brownian  motions,  we  know  from  Section  II,  that  the  generator  !£  is  given 


$£=  l  m  (x,t)  ~  +  \  l  b  (x,t)  - 

i-1  i  i, j=l  iJ  ^i^j 


(5.76(A)) 


where 


b.M  =  l  o..(x,t)o..  (x,t). 


ik  jk 


It  immediately  follows  that  the  generator  for 


dx  =  em(x,t)dt  +  /e  a(x,t)dB 


(5.77) 


becomes  eit ’,  where  is  defined  by  (5.76(A)). 
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Hence,  the  theorem  states  that  the  probability  distribution  generated  by 


where 


dy  =  em(y)dt  +  ✓  e  a  (y)dB, 


m(y)  =  lim  -  m(y,t)dt 

T+oo  Jq 


(5.78) 


—  1  r 
b  (y)  =  lim  - 


b  (y,t: 

J0 


(5.79) 


(j(y)  is  obtained  from  (b..(y)),  is  close  to  the  probability  distribution 
generated  by  (5.77)  on  the  internal  [0,T/e].  Examples  of  this  approach  for 
non-linear  systems  may  be  found  in  [5.37]. 

Perhaps  the  most  interesting  form  of  averaging  methods  can  be  found  in 
[5.36],  [5.38],  [5.39].  In  a  sense  these  constitute  an  extension  of  the 
central  limit  theorem,  for  diffusion  processes.  We  shall  mention  the  result 
in  the  form  stated  in  [5.36]  (which  puts  the  results  of  [5.38]  on  a  firm 
mathematical  foundation). 

Consider  the  standard  form, 


z(t)  =  eF(z(t) ,  x ( t ) ,  t),  z(0)  =  z 


where  z,F  are  n-vectors,  x  is  a  vector  stochastic  process,  where  for  each 
component  F^  of  F,  we  assume  that  there  is  a  constant  C,  such  that 

3F  92F 

i  F  .  |  <  C ,  |  - — -  |  <  C ,  |r - vH  <  C 

1  dZ  dZ  dZ 

j  j  k 


(5.80) 


(5.81) 


uniformly  in  (z,x,t). 


J 
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It  is  assumed  that  the  x-process  possesses  the  strong  mixing  property. 

This  essentially  means  that  as  time  intervals  [0,t],  [t  +  s,=°)  become 
further  separated  (i.e.  as  s  +  °°) ,  events  defined  on  the  x-process,  over 
these  two  intervals,  will  become  independent. 

Finally,  it  is  assumed  that  E{F(z,x(t) ,t) } =  0,  for  fixed  z,  and  the 
following  limits  exist  uniformly  in  z, 


lira  — 
T 


j  dtiJ  d'2  — 

n  n  J  x 


n  SF^z.xCt)  ,-t^) 


F  (z,x(t2)t2)  }  =  bi(z) 


lim  ~  dx  dx  Et F. (z,x(x  ) ,x  )F. (z,x(x  ) ,x  ) } 

T+oo  i-<0  J0  1  2  z 


Under  these  conditions,  the  process  z^  '( t )  defined  by 


Vz) 


(5.82) 


z^\t)  =  z(x/e2) ,  t  =  e2t. 


(5.83) 


converges  weakly  (i.e.  in  distribution)  to  a  diffusion  process  whose  generator, 
or  backward  operator,  j? ,  is  given  as 


-l  n  *\  2  ^ 

=  x  )  a  .  .  (z)  — - - —  +  \  b .  (z)  t-— 

2  •;  .,ij  3z , 3z .  .  ,  l  3z. 

i,J=l  i  l  J=1  j 


(5.84) 
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Since  we  can  always  associate  a  white  noise  coefficient,  Ito,  differential 
equation  with  the  generator  (5.84),  then  it  follows  that  the  statistics  of  the 
solution  process  of  (5.80),  will  converge  to  the  statistics  of  the  solution 
process  of 


dy  =  b(y)dt  +  G(y)dB, 


(5.85) 


where  the  vector  b(y)  and  matrix  G(y)  are  determined  by  the  (b^(z) ) , (a  (z) ) 
of  (5.82). 

One  important  note,  the  convergence  is  weak  convergence.  Therefore,  we 
cannot  say  that  the  actual  sample  response  behavior  of  (5.80)  will  approach 
the  response  of  (5.85),  we  can  only  say  that  their  statistical  properties  as 
governed  by  their  joint  probability  densities  will  become  close  as  e  approaches 
zero.  An  extension  and  further  development  of  these  ideas  are  due  to  Papanicolaou, 
[5.39]  . 

The  following  example  of  (5.39)  is  illustrative  of  the  procedure.  For 
further  applications  see  [5.40]  and  [5.19]. 

Consider  the  undamped  oscillator, 

x(t)  +  (tu  +  sn(t)  )x(t)  =  0,  (5.86) 


where  the  n-process  is  bounded,  with  covariance  y(x).  We  can  put 
(5.86)  into  the  usual  phase-space  form  as 
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Upon  applying  the  transformation 


(5.88) 


noting  that  the  matrix  exponential  e  will  generate  trigonometric  terms, 
the  equation  (5.87)  may  be  written  as 


Z1  2  z2  2 
en(t)  [—  sin  ut  +  — r—  sin  ut] 
2  u  2 

to 


2  2 

en(t)  [-z.  cos  ut  -  —  sin2utl, 
1  2u 


(5.89) 


Cne  final  transformation. 


z^  =  ercos0,  z2  =  -u)ersin6. 


(5.90) 


transforms  (5.89)  into 


r  =  en(t)G(0 , t) 


0  =  en(t)H(9,t) 


(5.91) 


where  G,  H  are  trigonometric  polynomials  in  (9,t).  Therefore,  G,H  are 
bounded  and  the  averaging  theorem  of  Khazminskii-Stratonovich-Papanicolaou 
can  be  applied  to  obtain  the  generator 


as;.-  ,  (**•»»» 
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3 


2 


+  b 


(a  +-|) 


where  from  (5.82)  one  obtains 


a 


b  =  wRef^l^-]  ,c  =  ulm  [^|^] 


(5.92) 


(5.93) 


with 


S(u) 


,  .  itur, 
y(x)e  dx. 

0 


Thus,  the  second  moment  approximation  to  E{x 
is 


2 

1 


(t) 


} 


on  the  internal 


[0,T/e2] 


O  *T  O  Cm  Cm 

Efx^(t)}  =  —  exp{^-  [ReS(2w)  -  2S(0) ] e  t}cos(2w-  —  ImS(2w))t 


+  exp[ojReS(2u)e2t] . 


(5.94) 


We  can  further  note  that  the  generator^  given  by  (5.92)  implies  that 
r,3  processes  are  independent.  (Sincej^?;  S£  +^Q) .  Therefore,  we  can  study 
the  r-process  through 


3r 


(5.95) 
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whose  associated  Ito  equation  is 


dr  =  bdt  ±  /b  dB. 


(5.96) 


Since  the  constant  term  b  is  positive,  it  follows  that 


lim  r (t)  =  °°  with  probability  one. 


(5.97) 


This  will  imply  that  (z, ,z»)  and  therefore  (x  ,x.)  are  growing  in  an 

12  2  1  2 
unstable  manner  on  the  internal  [0,T/e  ]. 

This  concludes  our  discussion  of  averaging  methods. 


SHE 
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VI.  Perturbation  methods* 


Consider,  for  example,  the  rth  natural  frequency  uj^  of  a  system. 

Assume  the  parameters  depend  upon  random  variables  X  ,  ...,  X  .  Then, 

1  m 

we  seek  an  expression  for  co^,  now  regarded  as  a  random  variable,  in  the 
form 


m  in  id 

u>  -  0)  +  z  +  Z  E  A  X,X4  +  ...  . 
r  i-1 1 1  i-u-1  1 J 

where  the  X^  are  regarded  as  small  perturbation  terms,  represents  the 
rth  natural  frequency  of  the  mean  system,  and  the  A^,  A  ,  ...  are  to  be 
determined.  While  some  authors  replace  X^  with  eY^,  e  being  the 
perturbation  parameter,  we  shall  not  usually  do  this.  Once  we  know  the 


we  can  obtain  statistical  properties  of  id 


r 


or  any  other 


quantity  of  interest.  Let  us  consider  a  general  formulation  of  this 
problem,  considering  natural  frequencies  and  normal  modes  first. 

Zarghame  [27]  sent  one  of  the  authors  a  method  of  this  class  which 
appears  well  suited  to  computation  and  which  makes  a  useful  suggestion 
on  how  to  introduce  random  parameters  that  merits  attention.  This 
method  has  never  appeared  in  print  in  so  far  as  we  know.  Therefore,  we 
shall  write  out  the  details  in  order  to  have  it  before  us. 

We  are  concerned  in  this  subsection  with  the  free  motion  of  a 
conservative  system.  Thus,  in  (6),  C  *  0,  f  *  0.  Now  the  elements  in 
the  symmetric  stiffness  matrix  K  are  determined  by  the  bars,  beams, 
columns,  joints,  etc.  making  up  the  structure.  The  uncertainties  in 

*  References  in  Sections  VI-IX  are  given  at  the  end  of  Section  IX. 
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the  structure  reside  in  these  elements.  Let  there  be  m  structural 

,th 


elements  in  the  structure,  and  let  the  stiffness  matrix  of  the  i 


structural  element  in  terms  of  . . .  be 

1  n 


\  «  (1  +  V^L  ‘ 


i  ”  1,  m  , 

which  produces  the  (nxn)  random  stiffness  matrix 


(6.1) 


<V  • 


(6.2) 


The  random  variables  X  ,  ...»  X  describe  the  uncertainty  present  in  the 

1  m 


structural  elements  and  we  assume 


EXi  =  0  ,  Var  85  . 


(6.3) 


tli  T 

K  is  the  mean  stiffness  matrix  of  the  i  element,  K  =  K  ,  i.e.  K  is 


symmetric  in  the  K^»  is  random  stiffness  element  corresponding 


to  q j  and  q^,  and  we  assume  masses  of  the  elements  do  not  change.  The 


advantage  of  (6.1)  and  (6.2)  is  that  statistical  dependence  of  the  K 


jk 


is  brought  in  a  straight-forward  manner.  We  note  also  that  we  can  write 
(6.2)  as 


K  =  K  +  EX.K 
—  i— 1 


IK. 


(6.4) 


which  gives  also 


EK  -  K  and 


9K 


9X. 


=  K. 


(6.5) 


where  K  is  the  stiffness  matrix  of  the  structure  with  each  number  taking 
its  mean  stiffness. 

We  can  now  write  (1.26)  as 


Iq  +  Kq  =  0  , 


(6.6) 


99  - 


where  q  is  the  (nxl)  column  vector  with  transpose 

T 

9  m  {q^ . 9n}  •  (6.7) 

Assume  normal  mode  motion 

q  ■  a  cos(wt  +  <f)  (6.8) 

with  a.  the  (nxl)  column  vector  defined  by 

aT  =  {a^ ,  ....  aj  .  (6.9) 

Then,  substituting  (6.8)  into  (6.6),  we  obtain 

(K  -  u>2I)a  =  0,  (6.10) 

where  again  I  is  the  (nxn)  unit  matrix. 

2 

The  squared  natural  frequencies  are  determined  by  the  n  roots  of 
the  equation 

det |K  -  u>2l|  -  0  ,  (6.11) 

2 

revealing  that  the  and  are  random  variables  since  K  contains 

random  variables.  Let  the  random  mode  corresponding  to  oi^  be  the  (nxl) 

column  vector  a  .  Then  we  can  write 
r 

(K  -  a)2I)a  -  0  ,  (6.12) 

r  r 

with  the  usual  orthogonality  relations 

aTla  =  0  ,  aTKa  »  0  if  s  *  r 

r  s  r  s 

aTIa  *  1  ,  aTKa  *  w2  ,  (6.13) 

rr  rrr 

where  "super  T"  denotes  transpose. 

We  are  now  interested  in  expressing  the  random  variables  o>^  and  a 
in  terms  of  a  series  in  the  random  variables  X^. 


The  expressions  we 


100  - 


seek  will  be  power  series  in  the  X^.  Notice  that  in  this  formulation  w^ 
and  are  regarded  as  functions  of  the  uncertainty  in  stiffness  in  each 
Individual  structural  member. 

Differentiate  (6.12)  with  respect  to  X^: 


3oi  „  3a 

-  2“r  Ixf  »“r  +  <K  *  V>  Ixf  '  0 


(6.14) 


Next  premultiply  (5.14)  by  a^  obtaining 


“r<ii  -  2“r^  I)ar  *  0  - 
T 

since  by  the  symmetry  of  K  (K  *  K  )  and  (5,12) 


aT(K  -  m2I)  -  0  . 
r  r 

Thus,  with  the  last  two  of  (5.13) 


3(0  ,  _ 

r  IT 

=*  ~  <»  K.a  . 
3X^  2u  r-i  r 

i  r 


(6.15) 


(6.16) 


This  is  to  be  evaluated  at  X,  =  ...  =  X  =  0  (i.e.  X  *=  0);  we  obtain 

1  n 


<!?>„  '  S'  2&A  . 

i  -r 


(6.17) 


where  the  underbarred  quantities  are  to  be  evaluated  for  the  system  with 
mean  stiffnesses.  We  note  that  (6.17)  give  the  sensitivity  coefficients 
[19,25]  of  with  respect  to  the  X^. 


The  a^,  k  *  1,  ...,  n  span  the  coordinate  space;  hence,  we  may 


write 


ES^a 

j  rl  j 


(6.18) 


We  substitute  (6.18)  into  (6.14): 
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do)  .  . 

(K,  -  2“  -*r-  1)°  +  (K  -  =  0  . 


ri  j 


(6.19) 


Now  premultiply  (6.19)  by  a£,  obtaining 


\<*i  '  2“r  ^  I)ar  +  *£<*  -  ‘  °  ' 


or,  for  k  *  r  and  with  the  use  of  (6.13) 


(to2  -  td2)0^)  «  -aTK  a  . 
v  k  r'ri  k-i  r 


Differentiating  the  next  to  last  of  (6.13)  with  respect  to  gives 


_  3a 

V  3X^  ■  0  * 


T  Cr) 

which  on  premultiplying  (6.18)  by  a^I  then  demonstrates  that  8^  '  =  0* 


Thus , 


T 

“kVi 


2  2 

0),  -Ul 

k  r 


k  *  r  , 


(6.20) 


and  hence  from  (6.18) 


3a  a^K.a 

r  k-i  r 

aX  Z  2  2  k  ’ 

i  j  0)  -o> 
r  k 


(6.21) 


where  the  prime  on  £  means  that  j  does  not  take  the  value  r.  When 


evaluated  at  X  *  0,  we  have 


T 

<*.K,a 


( _ =  r> 

' 3X  'o  .  2  2 

i  j  w  -a) 


(6.22) 


where  again  the  under  barred  quantities  are  evaluated  when  members  take 


S  .*  i  /  .* 


E 


on  their  mean  stiffnesses.  We  note  that  (6.22)  gives  the  sensitivity 
coefficients  [28]  of  the  mode  shapes  with  respect  to  the  X  . 


Next  differentiate  (6.14)  with  respect  to  X^: 


2 

3co  3w  3  w  3oi  3a 

“2(3X^  -53^  +  “r  3Xi3XJ)Iar  +  (^i"2aJr  "3X^  X)  3X 7 


(6.23) 


32a 


3oj  3a  „ 

+  (^-2"r  lx;  !)  HTL  -  ^"r1*  IXjXx; 


0  . 


Next  employ  (34)  for  3a  /3X,  and  3a  /3X,  in  (36)  and  premultiply  by  a1; 

ri  rj  r  J  J  r  * 

these  operations  yield 


3u  doi 
r  r 


,2 

3  u) 


r  v  T. 


-2(<< +  \  + 

i  j  i  j  k  J 

3  oi 


-2o>  r'8^)aTIa,  +  Z'B^aVa, 


r  3Xt  £  "ij  rk 


ri  r-j  k 


-2u>  ~  E'eJ^c^la  . 
r  3Xj  k  li  r  k 


Since  by  (6.13)  and  (6.20), 


reri)aII\  -  r6rVarIak  ■ 0  • 


and  by  (6.20) 


u‘K.0^  =  6^)(^-<*>J),  =  ^)(o,2-o.2)  , 


-j  k  ~kj  v~k  k' 


we  f ind 


.2 

3  o) 

r 

3X.3X, 


i  “j 


lr7//R(k)R(r)  fi(k)ft(rk,  2  2. 

2 ll  (6ri  ekj  +6rj  \i  )<VWr) 


3o>  3oi 

-2  — ^  ^1 
*•  5v  Tv  j 


3x1  ixy 


or  at  X  =  0, 


(6.24) 


(6.25) 


i 
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/  r  \  _  lf.  rn(k)n(r)  (k)  (r)Vn2  2. 

(3x75xT)o  2lZ/iri  6kj  +4j  0ki  (— k  — r 
1  J  k 


(6.26) 


3a»  3a> 

“^Ix^o^'SxjV 


To  go  further,  we  must  first  evaluate  3  a^/Sx^SX^.  Let 


£6°^  A. 


Ci3Xj  k  r,ij  k  * 


(6.27) 


since  as  with  (6.18)  the  a^  span  the  coordinate  space.  The  substitution 


of  (6.27)  into  (6.23)  yields 


3o)  3(0  3  0}  3(i>  3a 

'2(3x7  "3x7  +  “r  3x73X“)Iar  +  (V2“r  3xf  X)  ST.  (6*28) 

i  j  i  j  1  3 

3“  3a  o 

+  ST  »  1x7  +  jA  -  0  . 


-j  r  3X  '  3Xt  '  r  r ,  1  j  k 

T 

Premultiplication  of  this  equation  by  a^,  JL  *  r,  and  the  employment  of 
similar  relations  as  used  to  obtain  (6.25)  ultimately  yields  for  1^  *  r. 


8U)  ,  1  r  T  ^i  ^r 

8r,ij  2  2  ”l  -i  2  i  3X. *  3X. 

*  J  w  -w,  —  i  j 

r  1  J 


(6.29) 


3u>  3a 

+  W*"r  3X7>  3X71  • 

-  3  1 


If  we  differentiate  the  equation  before  (6.20)  with  respect  to  X^,  we 


S2a 

_T  _ r_ 

r  SX^X^ 


3aT  3a 


3Xj  3X1 


Then  on  premultiplying  (6.27)  by  a^i  and  employing  this  relationship  the 


result 


N  «  lVW.x  w  -  hVj.VWm''  -  ».  r  *  •  *  *  *  •  *>  .*■  *  - 
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(r) 


3a 


3a 


i.ij 


3X, 


3X, 


(6.30) 


J  i 

follows.  Equations  (6.29)  and  (6.30)  can  now  be  evaluated  at  X  =  0,  and 


we  obtain 


32a  (  . 

, _ r  V  _  rfl(k) 

\  IN  v  “  P  ,  .  • 

3X13XJ/o  k-r,ij-k 


(6.31) 


We  could  rewrite  (6.29)  and  (6.30)  employing  previously  determined 
expressions  for  the  first  partial  derivatives  contained  therein  but  this 
is  not  particularly  helpful.  The  procedure  for  going  further  is 
straight  forward  but  we  shall  not  write  out  the  details  in  order  to 
conserve  space.  Let  us  summarize  our  results  up  to  this  point. 

We  have  for  the  random  variable 


3o) 


a2 

3  o> 


r  *  -r  +  +  te'ixpifT’oVj  +  • 


-ij  -“i’“j 

where  the  partial  derivatives  are  supplied  by  (6.17)  and  (6.26). 


(6.32) 

We 


also  have  for  the  random  variable 


3a 


32a 


\  ■  5, +  ?<icVi  +  +  •••  • 


(6.33) 


i  i 

where  (6.22)  supplies  the  first  partial  derivative,  and  (6.29),  (6.30) 


supplies  the  derivatives  in  the  double  sum.  Before  looking  at  the 


statistics  of  m  and  a  ,  let  us  consider  how  the  K.  and  K  are  computed, 
r  r  — i 


Consider,  for  simplicity,  a  plane  frame  consisting  of  pin-ended 


,th 


straight  bars.  Let  the  10  bar  connect  joints  13  and  19,  for  example. 


Let  this  bar  have  length  1  ,  mean  area  A^,  mean  modulus  E^,  and  mean 


direction  cosines  and  as  shown.  Let  the  elastic  displacements 


I 
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at  the  joints  be  (u^,  v^)  and  (u^,  v^). 


The  mean  elastic  potential  energy  stored  In  the  bar  is 


„  1  A10E10  t2 
*  ”ilo  S'° 


(6.34) 


where  to  first  order  in  small  displacements 


610  =  (u19“U13)X10  +  (V19_V13)P10  * 


(6.35) 


The  substitution  of  (6.35)  into  (6.34)  gives,  on  arranging  terras, 


columns 


Qi 

u13 

V13 

•  •  • 

U19 

V19 

S  u13 

AEX2 

1 

AEXp 

1 

•  •  • 

AEX2 

.1 

AEXp 

1 

AEXp 

AEp2 

AEXp 

AEp2 

U  V13 

_1 

1 

e  •  • 

_1 

1 

• 

• 

• 

•  •  • 

• 

• 

^  rows 

• 

• 

•  •  • 

• 

• 

• 

• 

• 

•  •  • 

• 

• 

K 

AEX2 

AEXp 

AEX2 

AEXp 

w  u19 

1 

1 

•  •  • 

1 

1 

K 

_  AEXp 

AEp2 

AEXp 

AEp2 

b8  v19 

_1 

1_ 

•  •  • 

1 

1 

"*!<)• 


(6.36) 


where  all  other  entries  in  the  K  matrix  are  zero,  and,  hence 
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(1  +  xio*io  • 


(6.37) 


Thus,  according  to  (6.2) 


K  -  Z  K  . 
i-1  1 


(6.2) 


Note  that  X^  represents  the  combination  of  all  the  random  variables 

that  may  enter  E^,  1^q»  ^q*  and  W^q,  assum^n8  uncertainty  comes 

only  from  the  bar  and  not  the  joints.  If  the  major  share  of  the 

uncertainty  in  stiffness  of  bar  10  comes  from  the  joint  connections, 

then  XjQ  must  describe  this  fact;  in  this  case,  X^q  may  be  dependent  on 

the  random  variables  associated  with  those  bars  sharing  the  joints  with 

bar  10.  Finally,  it  may  be  advantageous  to  introduce  random  variables 

associated  with  only  joint  behavior  if  its  behavior  has  a  substantial 

influence  on  stiffness  of  the  structure.  We  conclude  from  this  brief 

discussion  that  the  X,,  ...,  X  may  be  independent  random  variables  if, 

1  m 

for  example,  only  bar  stiffness  need  be  considered  and  the  bars  do  not 

interact  with  each  other;  however,  it  is  possible  the  X^,  ...,  X^  may 

be  dependent  if  the  bars  interact  through  joint  behavior.  Equation 

(6.36)  demonstrates  that  by  attaching  random  variable  to  physical 

element's  stiffness  coefficient,  statistical  dependence  in  K  is  easily 

included  whether  or  not  the  X.,  ....  X  are  dependent.  Let  us  consider 

1  m 

the  statistics  of  etc  ... 

Consider  equation  (6.36)  first.  We  have,  on  taking  expectation, 


E,“r>  •  2r  +  l”'ax^>,EfXlV  +  -  • 


(6.38) 


Even  if  the  X's  are  independent  E{o>  }  *  u>  ,  since  the  EX,  *  0  terms  are 

r  -r  j 

still  present.  Now  square  (6.32)  and  take  expectation: 


8 


I 
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g 

B 

8 

R 


E{u)2}  -  to2  +  2o)  EE(- 
r  -r 


.2 

9  (i) 


3Xi3Xj)oE{XiXj} 


(6.39) 


3(0  3(0 

+  “<^>„<^>„I,xiV 

3(0  3  (o 

32  32 

+  +  •••  • 


We  can  now  approximate  Var  it  is  defined  as 


Var  u>2  =  E{uj2}  -  [E{o>r> ] 2  .  (6.40) 

Thus,  it  is  a  straight  forward  task  to  approximate  the  first  two  moments 

of  go  . 
r 

If  we  extend  (6.32)  to  cubic,  quartic,  ...  terms  in  the  X^,  then 

(6.38)  and  (6.39)  would  contain  additional  terms.  How  far  we  should 

continue  this  process  will  depend  on  relative  size  of  the  terms 

containing  EtX^X^},  EtX^X^X^},  etc.  and  what  information  we  have  that 

would  enable  as  to  evaluate  these  expectations.  It  is  not  usual  that  we 

can  evaluate  any  more  than  E{X_^Xj}. 

Let  us  consider  the  determination  of  the  distribution  of  <o  .  To  do 

r 

this,  we  require  the  joint  distribution  of  X  ,  ...,  X  ;  denote  the  joint 

1  m 

probability  distribution  by  f  (x. . x  ).  Then, 

mi  m 


F  (w)  =  P{u  <  u>}  (6.41) 

r  r 

*  J • • • f f  (X|,  ...,  x  )dx,  ...  dx 
mi  ml  m 

where  the  multiple  integral  is  over  all  x^  such  that  <  u >,  and  where 
from  (6.32), 


-  108  - 


do>  .  3  w 

“r(Xl . V  Sr  +  +  S'iSTTWj  +  -  (6-42> 

i  i  ij  i  j  J 

Evaluation  of  (6.41)  is  a  formidable  task  even  for  m  reasonably  small. 
This  paragraph  is  thus  largely  cultural  in  so  far  as  practical 
application  is  concerned. 

We  have  put  in  this  detailed  treatment  of  Zarghame's  method  since 

it  is  not  in  the  literature  and  appears  useful  as  mentioned  earlier.  In 

particular,  we  note  that  it  gives  sensitivity  coefficients  for  natural 

frequencies  and  correspondence  normal  modes  plus  series  expansions  for 

these  quantities  in  terms  of  the  random  variables  X  ,  ...,  X  that 

1  m 

define  the  uncertainty  present  in  the  stiffness  matrix  K.  Moments  of 


the  quantities  are  easily  obtained,  but  it  is  practically  impossible  to 
obtain  distributions  for  the  natural  frequencies  and  corresponding 
normal  modes.  For  confidence  interval  location  and  size  for  a  natural 
frequency,  for  example,  we  must  approximate  using 


E{u  }  ±  3,  Var  oj 
r  \  r 


(6.43) 


as  a  rough  indication  of  a  99%  confidence  interval.  This  interval  gives 


us  some  idea  of  the  spread  in  a  natural  frequency  and  it  could  be 
employed  to  make  reasonably  sure  that  no  steady  excitation  frequencies 


were  contained  therein  for  all  cu^.  Alternatively,  we  might  employ  the 
signal  to  noise  ratio 


E{ur> 

|Var  id 


(6.44) 


to  obtain  an  idea  of  how  important  stiffness  uncertainty  is  for  natural 
frequency;  if  (6.44)  is  greater  than  20  or  30  say,  we  would  regard  the 
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location  of  oi^  as  deterministic;  on  the  other  hand,  if  (6.44)  is  less 
than  5-10,  it  might  be  unwise  to  ignore  this  level  of  variability  in  the 
location  of  depending,  of  course,  on  the  consequences  of  such 
uncertainty. 

Next,  let  us  consider  frequency  response.  The  Fourier  transform  of 
(1.26),  with  the  unit  mass  matrix  I  replaced  by  M,  is 

(K  -  oi2M  +  iu>C)Q  -  F  ,  (6.45) 

where 


Q  -  2 7  /  qeia,td“*  F  -  W  /  fei"td‘ 


2ir 


(6.46) 


Let  the  (nxn)  matrix 


(K  -  or  M  +  iwC)  *  Z(o>)  «  Z  . 
The  inverse  of  Z,  written  Z  1  ( to)  »  Z  *,  satisfies 


(6.47) 


ZZ  1  -  z-1z  -  I 


(6.48) 


Writing  (6.45)  as 


ZQ  -  F  , 


(6.49) 


we  find  on  using  (6.48)  that 


i 

I 

3 

I 

1 

I 


q  -  Z-1F  ,  (6.50) 

where  the  explicit  dependence  of  Q,  F,  and  Z  *  on  to  is  not  shown. 

The  matrix  Z  *  is  called  the  frequency  response  matrix  of  the 
system.  The  physical  meaning  of  Z  *  is  as  follows:  write 


Z-1(u)  -  {Z“^(t0>>  ; 


,-l 


then  the  element  Z  is  the  complex  response  amplitude  at  q  due  to  a 
3*  J 


no  - 


unit  force  (with  proper  dimensions) 


iut 


acting  at  qfc;  i.e. 


v  iwt 

V  >* 

is  the  response  at  q  (output)  due  to  the  above  force  at  q  (input). 

J  k 


-1 

jk 


The  symmetry  of  the  matrices  M,  K,  and  C  establishes  the  symmetry  Z 
-1  -1  T  -1 

Z  ..or  (Z  )  ■  Z  .  In  a  lightly  damped  structural  system,  when 

^  J 

amplitude  |Z^(ai)|  is  plotted  as  ordinate  against  w  as  absicssa,  sharp 

peaks  will  appear  in  these  curves  in  the  neighborhood  of  natural 

frequencies  of  the  undamped  system.  This  means  that  resonance  (large 

amplitudes)  occurs  in  at  least  some  of  the  q  due  to  this  force  at  q., 

J  k 

when  a)  is  near  to  one  of  the  undamped  natural  frequencies  of  the  system. 
Put  another  way,  if  we  regard  the  system  as  a  mechanical  filter,  only 
frequencies  close  to  the  undamped  natural  frequencies  where  there  is  a 
peak  in  |Z  ^(m)|  will  show  up  in  the  output  q  for  the  input  q.  . 

Knowing  Z  (w)  we  obtain  the  response  vector  q  by  taking  the 
inverse  Fourier  transform  of  (6.50);  thus, 


q  -  /  eiu,t  Z_1(u>)F(u))du>  .  (6.51) 

If  the  excitation  vector  f  is  a  wide  sense  stationary  random  process,  or 
if  f  is  periodic,  Z  *(o>)  is  also  the  quantity  needed  to  obtain  the 
response  q.  Hence,  our  interest  in  Z  * . 

The  inverse  Fourier  transform  of  Z  ^(u>)  produces  the  impulse 
response  matrix  H(t),  a  (nxn)  matrix: 
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H(t)  -  /  eiWtZ  ^ (o))dto  ,  t  >  0  , 


(6.52) 


t  <  0  . 


The  element  h  (t)  Is  the  response  at  q  when  a  unit  velocity  is 

J 


produced  only  at  q^,  for  q^  *  ...  *  q^  °  0  at  t  *  0.  Again  by  symmetry. 


h  ^  =  h^.  By  the  convolution  theorem  of  Fourier  transforms,  we  may 


replace  (6.51)  by 


q  -  /  H(t-T)f(r)dx  . 


(6.53) 


Equation  (6.52)  shows  that  the  impulse  response  matrix  H(t),  which  is  in 


the  time  domain,  is  equivalent  to  the  frequency  response  matrix  Z  (u>), 


which  is  in  the  frequency  domain. 


We  noticed,  while  writing  out  Zarghame's  expansion  method  for 


natural  frequencies  and  normal  modes,  that  it  is  possible  to  extend  it 


to  include  frequency  response. 


We  assume  Z  is  given  by  (6.47): 


Z  -  K  -  w  M  +  iu>C  . 


(6.47') 


Let  X  ...,  X  be  the  set  of  random  variables  that  describes  the 
1  tn 


distributions  of  member  properties  as  contained  in  K,  M,  and  C.  Now  m 


is  the  total  number  of  distinct  sources  of  variability  in  K,  M,  and  C, 


which  may  be  bigger  than  the  number  of  members  in  this  case.  (For 


example,  parallel  acting  but  distinct  sources  of  damping  and  stiffness 


in  the  same  number  would  require  two  different  X's.)  We  recall  that 


before  m  equals  at  most  the  number  of  members  according  to  Zarghame's 


formulation. 


We  seek  an  expression  of  Che  type 


,-l 


z"1  +  z  (3Z 


-1 


J 


i  aV1 
3Xj~^oXj  +  IT  “ 


7\)oxA+  -  • 


(6. 54) 


We  also  seek  a  simple  method  for  evaluation  of  the  partial  derivatives. 
Consider 


ZZ_1  -  I  . 


(6.48') 


The  differentiation  of  this  equation  with  respect  to  X^  gives 


3Z  „-l 


az 


-l 


ax , 


0  ; 


(6.55) 


-1 


premultiplication  by  Z  produces 


3Z 

1* 


-1 


_7-i  az  -i 
z  "5x"~  z 


(6.56) 


i  i 

Since  by  (6.47'),  noting  that  w  is  a  parameter. 


*  K,  -  w2M,  +  iuC  , 

Xj  _j  _j 


(6.57) 


we  now  find 


(||— )Q  “  +  1^)Z_1  ,  (6.58) 

where  under  bar  means  evaluation  is  for  the  system  with  all  member 
parameters  taking  their  mean  values.  Equation  (6.58)  is  the  sensitivity 
coefficient  of  the  matrix  Z 


Next,  the  differentiation  of  (6.55)  with  respect  to  X^  produces 


a2z  „-i  ,  az  az"1  ,  az  az"1  ,  „  a2z_1 

Again,  premultiply  by  Z  ^  and  rearrange  to  obtain 


0 


(6.59) 
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is 


32Z~‘  ,-l  32Z  -1  3Z  iz'1  ,-l  3Z  3z" 

We  now  use  (6.56)  and  (6.57)  to  yield 


(6.60) 


32Z_1 


.-1. 


,-l. 


<1x^0  (Kj-  Mj+icoC^Z  (^-o,  M^-C^Z 

+  Z  1  (K^-(U2M^+ia)C^  )Z~ 1  (  K  j-a)2M  j+iwC  ^  )Z.' 
since  also  by  (6.57) 


-1 


(6.61) 


,-l 


a2z 


dX 


A 


=  0 


(6.62) 


It  is  straight  forward  to  obtain  the  higher  derivatives  of  Z~*.  Thus, 
(6.54)  can  be  carried  as  far  as  needed. 

The  first  obvious  advantages  of  this  formulation  is  that  it  applies 
to  any  type  of  structure  with  variability  in  K,  M  and  C.  Second  the 
needed  information  is  contained  in  the  system  with  mean  elements  only; 
once  this  information  is  in  hand,  everything  else  follows.  Because  of 
(6.62),  the  next  differentiation  of  (6.59)  with  respect  to  X^,  say,  will 
have  no  second  or  higher  derivatives  of  Z,  and  hence,  to  proceed  further 
is  not  difficult.  It  follows  that  this  formulation  appears  attractive. 
Convergence  requires  attention,  of  course.  Further,  numerical 
computational  ease  must  be  established,  case  in  [44]  although  the  ideas 
there  expressed  are  interesting.  Normal  coordinates? 

Randomness  in  parameters  can  be  introduced  into  the  coefficients  in 
K,  M,  and  C,  or  through  structural  elements  as  suggested  above,  or  by 
assuming  that  the  natural  frequencies  are  themselves  random  variables  as 


in  [35,  36]. 


The  purpose  of  this  short  digression  is  to  emphasize  to  the  reader 


that  the  choice  of  the  equations  of  motion  merits  careful  consideration 
so  that  what  is  important  in  a  problem  is  included.  Let  us  now  return 
to  the  above  references  [1,9,11,21,24,27,28,30,31]. 


Lord  Rayleigh  [1]  was  not  interested  in  our  topic.  However,  he 
did,  for  a  different  reason,  find  formulas  for  the  deterministic  change 
in  a  natural  frequency  and  corresponding  normal  mode  when  the 
coefficients  in  T  and  V  (in  the  normal  coordinate  formulation)  are 
subject  to  small  deterministic  changes.  His  formula  for  the  changed 
natural  frequency  is  (in  our  notation) 


_  k  +5k  n  (5k  -o)^6m  )^ 
j2  =  — r  rr  _  rs  -s  rs 

r  m  +5m  _ ,  ,2  2, 


U  i  /  L.  i.  v 

rr  s-l  m^mr(  0)^-0)^) 

where  s  *  r  in  E'.  Randomness  can  be  introduced  by  replacing  5k  and 


rs 


^mrs  smaH  random  variables.  The  first  term  represents  the  change  in 


u>r  due  to  change  in  mass  and  stiffness  without  changing  the  mode  shape; 


the  second  term  is  due  to  the  change  in  mode  shape.  This  formula  is  not 
used  today,  since  changes  in  parameters  in  the  normal  coordinate 
formulation  are  not  of  direct  concern.  We  bring  it  in  because  it  shows 


that  this  master  of  small  iteration  was  aware  that  small  parameter 
,2 


changes  may  alter  by  large  amounts,  which  is  of  concern  today. 

Reference  [9]  is  concerned  with  the  sensitivity  coefficients  of  the 
buckling  load  of  plates  with  random  thickness  and  []  discusses  the 
vibration  and  buckling  of  a  column  with  spring  supports  and  axial  loads 
treated  as  random  variables  and  with  material  and  geometric  properties 
considered  as  correlated  random  functions.  The  key  formula  in  [9], 
attributed  by  the  author  to  Jacobi,  is  the  same  as  (36.17),  which  is  the 


interesting  point.  Reference  [31]  employs  a  perturbation  method  (only 
linear  term  in  (6.32))  for  eigenvalue  change;  the  interesting  point  is 
that  a  Monte  Carlo  simulation  is  employed  to  check  the  analysis.  The 
computer  simulation  is  briefly  described.  The  results  from  analysis  and 
simulation  diverge  as  the  variance  in  the  parameters  increase  and  this 
is  disp^  ved  graphically. 

References  [21,24,28,31]  employ  either  sensitivity  coefficients  or 
linear  term  perturbation  expansions  (first  two  terms  on  right  of  (6.32) 
and  (6.33))  to  examine  influence  of  random  parameters  on  natural 
frequencies  and  normal  modes.  Their  techniques  are  in  the  same  general 
form  as  given  above  following  [27].  Except  for  the  last,  randomness 
resides  in  the  terms  in  K,  M,  C;  in  the  last,  the  component  mode 
synthesis  method  [16]  is  used  directly  to  derive  the  perturbation 
equations,  which  makes  this  paper  [31]  potentially  interesting  to  those 
confronted  with  an  actual  problem,  and  hence  randomness  may  reside  in 
the  structural  elements.  Numerical  difficulties  in  carrying  out  the 
computations  are  discussed  in  [28]  and  [30].  Reference  [31]  mentions  a 
computer  code  of  NASA. 

The  linear  chain  geometry  assumed  for  the  physical  system  in  [11] 
makes  it  possible  to  employ  a  different  technique  to  derive  the 
perturbation  expansion  than  described  above.  This  technique  employs  a 
transfer  matrix  method  [68,69]  and  it  is  applicable  whenever  a 

The  references  [1,9,11,21,24,27,28,30,31]  address  the  eigenvalue 
(natural  frequency)  and  eigenvector  (normal  mode)  problem  in  structural 
systems  by  perturbation  methods.  Before  discussing  methods  or 
techniques  involved,  it  is  important  to  understand  at  the  outset  that 


the  geometry  of  the  structure,  how  its  equations  of  motion  are 
assembled,  the  final  mathematical  form  of  the  equations  of  motion,  and 
how  randomness  in  parameters  is  introduced  have  a  profound  influence  on 
the  nature  of  the  results  obtained. 

A  structure's  geometry  can  be  in  the  form  of  a  linear  array  (chain) 
of  elements  that  may,  for  example,  consist  of  simple  harmonic 
oscillators  strung  together  in  a  line,  beam  segments  continuously 
connected  at  a  sequence  of  supports  in  a  line,  etc.  The  geometry  is  the 
simplest  possible  in  such  arrangements.  Plate  or  shell  type  structures 
have  a  two  dimensional  grid-like  geometry  and  are  next  in  order  of 
complexity.  Finally,  we  have  the  general  case  in  which  one  or  two 
dimensional  geometries  are  either  missing  or  are  interconnected  in  a 
complex  manner. 

The  equations  of  motion  depend  on  the  coordinate  choice, 
particularly  when  the  fact  that  mass  is  always  distributed  is  taken  into 
account.  Reference  [26]  discusses  methods  of  making  this  choice  and 
illustrates  the  substantial  difference  in  response  that  can  occur  due  to 
different  choices.  Reference  [16]  also  discusses  a  component  mode 
synthesis  method  for  selecting  coordinates  and  assembling  the  equations 
of  motion.  We  cannot  present  any  of  this  material  here  in  spite  of  its 
importance . 

A  coordinate  transformation  of  the  equations  of  motion  is  sometimes 
employed  as  in  [9,  44].  The  altered  form  of  the  equations  may  be 
advantageous  for  our  purpose  as  in  [9]  bu  this  does  not  appear  to  be 
the  structure  consists  of  a  chain  of  cells  or  units.  Let  us  consider 
this  method  in  some  detail. 
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Consider  the  system  shown  in  Figure  which  is  linear  chain  of 


oscillators  (See  [11])  and  where 
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we  have,  for  the  moment,  ignored  randomness  in  the  parameters.  The 


kinetic  and  potential  energies  are 


2T-TY3  ' 


(6.63) 


2V  «  t  k/Uj-Uj.;)  . 


which  gives  as  equations  of  motion 


mjU  j  +  kj^Uj_Uj-l^  ~  kj+l^Uj+l”Uj^  =  J  •••>  n_1  (6.64) 


mu  +  k  u  =  f  e  . 
n  n  n  n  o 


Introduce  a  new  coordinate 


wj  =  Vi  ■  uj  ; 


(6.65) 


iwu  r. 

u  =  x.e  ,  x  =  0  , 

j  3  o 


(6.66) 
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and  let 


We  now  write  the  equations  of  motion  as 
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where  I  is  the  2x2  unit  matrix 
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(6.68) 


(6.69a) 


(6.69b) 


We  can  now  relate  the  first  displacement  vector  d  to  the  last  d  by 

o  n 
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d  -  n  (I+T  )d  . 


(6.70) 


The  matrix  I+T  transfers  d  to  d  ;  hence,  the  name  "transfer 
J  J+1  J 

matrix".  Any  structure  whose  geometry  is  a  linear  array  of  units,  such 

as  in  the  above  Figure,  or  of  beams,  etc,  can  be  treated  in  this  manner. 

n 

Let  the  elements  in  the  (2x2)  matrix  IT ( I+-T  )  be  a.,  (to).  Then 

J 

(6.70)  becomes 


aUW  a12(a°\/V 

a2.(u>)  a22(u>)l  fo 

'  kn  i 


(6.71) 


The  natural  frequency  equation  of  the  chain  is  obtained  by  letting 
f  =0  and  taking  the  first  equation  in  (6.71);  it  is 


au(»)*n  -  0 

Since  x^  cannot  be  zero,  the  frequency  equation  is 


(6.72) 


ail(u>)  «  0  . 


(6.73) 


In  [11],  details,  which  do  not  concern  us  here,  are  worked  out  for  the 
natural  frequencies  when  k^  =  k,  =  m. 

Now  assume  the  masses  are  random  variables  taking  the  form 


nij  *  (1  +  X^)m  ,  (6.74) 

where  is  a  dimensionless  random  variable  with  mean  zero.  Introduce 


+  o 
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then,  with  k^  ■  k, 


(6.75) 
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I  +  Tj“I  +  T  +  Ej  , 


(6.76) 


Now  (6.70)  takes  the  form 


d  -  II  (I+T+E  )d 
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n 


(91) 


or 


d  -  [(I+T)n  +  Z  (I+T)j~1E.(I+T)n“J 
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(6.78) 


+  E  E  (I+T)J_1E  (I+T)k_;,_1EL  (I+T)n_k 
j-lk-j+1  J 
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A  number  of  substitutions  then  make  it  possible  to  expand  a >^_  in  a  series 


in  the  X^: 


n 


n-1  n 


u  ■  u  +  I  u.  4X.  +  E  »,.X*  +  E  E  a>  X.X,  +  ...  .  (6.79) 

r  “r  j-1  ^  *  j-1  2J  J  J-lk-j+1  ljk  A 


The  main  point  to  note  is  that  normal  modes  of  the  mean  system  are  not 
employed.  Let  us  contrast  this  method  with  that  of  Zarghame's. 


Return  to  (6.63),  let  k^  -  k  and  -  m  and  employ  (6.74)  to  obtain 


M  «  M  +  EMjXj 


K  -  K 


(6.80 


of  the  system  with  mean  elements  (i.e.  M  =  M,  K  *  K) .  There,  from  (A7) 


(6.82 J 
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We  observe  from  (6.82)  that  the  j  sensitivity  coefficient  of  oi 


depends  only  upon  w^,  and  M^.  Let 


T  ,  (1)  (2) 

o  *  (a  ,  a 
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)  . 


Then  (6.82)  becomes 


8V  *0*  ,  ( j),2 

(Sx^o - 2~  <2r  >  ; 


(6.85) 


remembering  from  (A6)  that  the  a'^  have  the  dimension  l/\|m,  we  see 
that  the  right  of  (6.85)  has  the  correct  dimension.  To  see  that  the 
sign  is  correct,  observe  from  the  first  two  terms  on  the  right  of  (6.32) 
that 


mm  .  .  . 

«r  "  ^r  "  E  (V y  Xj 


(6.86) 


Positive  Xj  mean  increase  in  mass.  We  know  from  Rayleigh's  Principle 


[1,  Sect.  ]  that  an  increase  in  mass  lowers  or  leaves  unchanged  every 
natural  frequency;  thus,  the  negative  sign  on  the  right  of  (6.85)  is 


correct. 
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The  importance  of  the  sensitivity  coefficients  (6.85)  and  (7A) 


resides  in  the  fact  that  they  reveal  by  their  magnitudes  those  u>^  that 


are  either  sensitive  or  insensitive  to  uncertainty  in  member  values. 


,00 


On  evaluating  the  0^  '  from  (6.83),  we  can  write  down  the  right  of 


(6.84).  The  main  point  to  observe  is  that  now  we  need  all  a  and  to. 

- r  — r 

Thus,  we  can  obtain  the  partial  derivative  values  appearing  in  (6.32). 

We  refrain  from  these  obvious  details  to  conserve  space. 

Let  us  contrast  Zarghame's  procedure  with  that  given  in  [11]. 

In  the  former,  we  need  all  and  to  to  obtain  (6.32);  in  the 
latter  we  need  all  to  but  not  the  normal  modes  since  the  intermediate 
coordinates  were  eliminated  in  obtaining  (6.70).  Thus,  Zarghame's 
procedure  does  require  more  information  than  required  in  [11].  However, 
since  all  computer  codes  produce  the  along  with  the  ui^,  the 
additional  information  required  in  the  former  is  always  available 
anyway.  Hence,  relative  to  effort  the  two  procedures  do  not  differ 
substantially. 


The  sensitivity  coefficient  in  (6.79)  is  given  by 


—  6.  .  cos  6 
m  1  j  o 


(6.87) 


where 


e. 


o  /_  _ \  =  (2r-l)r 

V  ’  5  2(2n-l)  » 

2  .  2 


(6.88) 


~  o—rr  sin  2 j  8  tan  6 
lj  2n+l  J  o  o 

While  these  are  easy  to  evaluate,  the  physical  significance  of  changes 
in  particular  parameter  values  Is  not  as  easy  to  grasp  in  (6.87)  as  is 
the  case  with  the  corresponding  formula  in  (7A).  This  difficulty 


increases  with  ^j*  an<*  ^ijn  Hence,  from  the  point  of  view  of 
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physical  understanding  of  what  is  of  significance  in  a  system  relative 
to  variability  in  natural  frequencies,  Zarghame's  procedure  appears 
superior  to  the  procedure  given  in  [11]. 

The  main  advantage  of  the  method  given  in  [9]  is  that  by  exploiting 
the  system  geometry  explicit  formulas  can  be  written  out  for  the 
quantities  of  interest.  However,  since  all  computations  are  now 
performed  on  a  computer,  this  is  no  longer  an  advantage. 

Two  other  points  deserve  comment.  First,  the  formula  for  in 
[30]  only  employs  the  first  two  terms  on  the  right  of  (6.32).  This 
makes 


E[(i>  }  *  ai  , 
r  — r 


whereas  in  (6.32) 


"V  ’  2r  +  1  “(^JS^WYk1  +  -  • 


Thus,  in  [30],  the  mean  of  is  the  mean  system  However,  from  the 

formula  above  this  is  not  correct  not  the  case.  Also,  (6.32)  gives  a 
different  formula  for  variance  of  than  given  in  [30].  It  follows 
that  the  method  given  in  [30]  is  incomplete.  Second,  sensitivity 
coefficients  also  play  an  important  role  in  other  types  of  system 
behavior  analysis  such  as  automatic  control  [14,15,18,19,25].  Since 
large  space  structures  will  contain  control  systems,  uncertainty  in 
control  system  parameters  coupled  with  uncertainties  in  the  structural 
system  parameters  must  be  kept  in  mind. 

References  [5, 8, 11, 20, 29, 35, 36, 37, 39, 40, 42, 43, 48, 49, 50, 52]  discuss 
frequency  response,  impulse  function  (Green's  function,  impulsive 
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admittance,  lmpedence)  mainly  by  perturbation  methods.  Those  that  do 
not  involve  perturbation  methods  are  [5,8,20].  Let  us  take  up  the 
latter  group  first. 

The  technique  involved  in  [5,8]  for  assessing  the  variability  of 
frequency  response  (or  impedance)  is  what  might  be  called  the  direct 
method.  For  a  one-degree  of  freedom  system,  we  have 

Z_1(a>)  -  - J -  .  (6-89) 

K-<*>  M+-iu»C 

where  K,  M,  C  are  scalar  random  variables.  There  are  cases  where 
knowing  the  distributions  of  independent  K,  M,  and  C  explicit  results 
can  be  obtained.  While  interesting  results  can  be  obtained  by  this 
method,  it  is  clear  this  technique  is  of  limited  practical  use  in 
complex  systems. 

Reference  [20]  starts  by  making  assumptions  about  the  statistical 
characteristics  of  the  Green's  function  G  of  a  system,  and  writes  the 
response  as 


q(t)  -  /  G(t,T)f(T)dT  . 
.00 

G  is  then  related  to  the  equation  of  motion 


by  assuming 


q  ■ 


-Aq 


(6.90) 


(6.91) 


A  -  A  +  N(t)  ,  (6.92) 

where  A  is  a  constant  matrix  and  N(t)  is  a  normal  noise  process.  On 
letting 


(6.93) 


t 

W(t,T)  -  /  N(0)d6  , 

T 

we  find 

G(t,T)  -  [A+N(T)]e"-(t-T)+W(t,T)  . 

Means  and  moments  of  q(t)  can  now  be  obtained  In  terms  of  the  moments  of 
G(t,x).  If  we  take  N(t)  to  be  independent  of  time,  then  this  technique 
applies  to  our  problem. 

We  observe  from  (6.94)  that  even  the  expectation  of  G  is  not  easy 
to  evaluate  because  of  N(t)  in  the  coefficient  and  W(t ,t)  in  the 
exponent.  By  assuming  N(t)  is  generated  by  filtering  a  white  noise 
process,  it  is  shown  that  the  moments  of  G(t,T)  can  be  evaluated. 

Further  investigation  of  this  interesting  technique  is  necessary  before 
we  can  determine  if  it  can  be  east  in  a  form  useful  to  us.  Even  if  it 
can,  additional  work  is  necessary  to  establish  that  it  can  be  applied 
when  N(t)  Is  constant  or  only  slowly  varying. 

Advantage  Is  taken  of  a  specific  structural  shape  from  the  outset 
in  [11,39,50].  The  first  two  assume  a  linear  chain  of  similar  elements 
differing  in  the  random  nature  of  element  parameters;  the  structure  in 
[50]  is  a  circular  chain.  Linear  one  degree  of  freedom  damped 
oscillator  elements  are  assumed  in  the  first  reference;  damped 
Bernoulli-Euler  beam  elements  with  random  lengths  are  assumed  in  the 
second;  and  a  continuous  distribution  of  linear  spring  connected  linar 
one  degree  of  freedom  oscillator  elements,  as  in  buckets  or  a  turbine 
disc,  is  assumed  in  the  third.  When  damping  is  assumed,  it  is  taken  as 
small  so  that  undamped  normal  modes  can  be  employed.  The  transfer 
matrix  technique  is  employed  in  [11,39]  since  a  chain-like  structure  is 


rvv' 
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assumed;  this  makes  It  possible  to  obtain  analytical  results  that  are 
complex  indeed.  Because  of  this  complexity,  it  will  require 
considerable  effort  to  obtain  numerical  results,  and  it  is  by  no  means 
clear  that  what  is  proposed  is  as  easy  to  use  as  (6.54).  The  periodic 
circular  structure  in  [50]  makes  it  possible  to  employ  Fourier  series; 
results  on  natural  frequencies,  normal  modes,  and  response  are  obtained; 
on  those  rare  occasions  when  a  structure  has  a  circular  form,  the 
technique  employed  could  be  considered  but  not  otherwise.  Excellent 
graphical  results  that  assist  in  understanding  in  a  qualitative  sense 
the  effect  of  disorder  on  response  are  presented  in  [11,39].  All  point 
out  that  high  variability  in  response  will  occur  in  lightly  damped 
disordered  systems.  The  techniques  employed  lack  the  generality  of 
(6.54)  and  are  not  of  direct  interest. 

Let  us  consider  [35,36,37,42]  next.  In  each  of  these  references 
random  parameters  are  introduced  in  special  ways  which  renders  the 
technqiues  used  and  results  obtained  of  limited  use.  [42]  does 
introduce  a  new  quantity  of  some  interest;  they  consider  the  equation 


mx  +  k( l+e)x  *  f (t) , 

where  e  is  random  with  zero  mean,  and  introduce 


V[x4(t)] 


E{x2(t)}-x2(t)( 

<x2(t)> 


(6.95) 


as  a  measure  of  the  time  for  the  mean  square  response  E{x4(t)}  to 

2 

deviate  from  the  unperturbed  response,  i.e.  x  (t)  with  t  =  0.  The 

2 

normalization  with  respect  to  the  time  average  <x  (t)>t  is  selected  so 
2 

that  V[x  ( t ) ]  ♦  1  as  t  ♦  *,  A  simple  expression  is  then  produced  for 
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the  envelope  of  V,  namely 

1  2  2  2 

2  -  T  U  °  C 

V[X  (t)3env.  "  1  '  e  »  <6’96> 

2  2 

where  u>  *  k/m,  o  *  Var  e.  V  is  essentially  the  growth  in  uncertainty 

2 

in  the  response  as  a  function  of  a  and  t.  This  is  a  nice  idea  that 
merits  development  since  it  says  that  when 
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response  location  is  lost. 

Papers  [48,49]  merit  attention  not  for  techniques  involved  but 
rather  for  some  qualitative  results  that  may  be  of  interest  in 
connection  with  large  space  structures.  Reference  [48]  discusses  wave 
propagation  in  long  beams  with  many  supports,  where  there  is  random 
variation  in  length  among  the  beam  segments.  The  point  of  interest  is 
that  the  random  variation  in  length  among  the  beam  segments  has  a 
substantial  influence  on  which  type  of  flexural  waves  will  propagate  arid 
which  will  attenuate.  Reference  [49]  is  concerned  with  the  confinement 
of  vibration  to  certain  regions  of  a  structure  due  to  structural 
disorder.  The  reason  for  noticing  them  resides  in  the  fact  that  in  a 
large  space  structure  it  may  be  desirable  to  introduce  structural 
irregularity  in  order  to  prevent  wave  motion  from  propagating  throughout 
the  structure  and/or  confining  vibration  to  a  favorable  region  of  the 
structure.  These  references  would  then  form  a  useful  starting  point. 

It  should  be  noted  that  the  "receptance  method"  [50]  ,  frequency  response 
method,  and  the  mobility  method  provide  essentially  the  same  approach; 
the  advantage  of  the  first  and  third  of  these  methods  lies  in  techniques 


for  obtaining  the  frequency  response  in  a  sequential  manner.  Referene 
[52]  investigates  failure  probability  in  a  structure  with  uncertain 
properties;  it  emphasizes  the  Importance  of  considering  these 
uncertainties  when  estimating  such  probabilities. 
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VII .  Llouville  Equation 

The  techniques  discussed  in  Section  5  are  based  directly  upon  the 
equations  of  motion.  For  example,  the  perturbation  expansion  of  the 
frequency  response  Z  *(w)  employs  the  equations  of  motion  to  derive 
Z  ^ (co )  and  the  relations  it  satisfies.  There  is  another  approach  based 
upon  the  Liouville  equation  for  the  time  evolution  of  the  joint 
probability  distribution  function  of  the  state  space  (2nxl)  column 
vector  x  and  the  system  parameters.  We  consider  that  approach  in  this 
Section. 

The  use  of  the  Liouville  equation  in  mechanics  and  statistical 
mechanics  is  of  long  standing  [51,52,53,  for  example  and  going  back  to 
Maxwell] .  These  references  do  not  assume  system  parameters  are  random 
variables,  and  average  of  quantities  under  equilibrium  conditions  is  of 
main  interest.  While  not  of  direct  interest  to  us,  it  is  possible  to 
adapt  these  early  methods  to  our  needs.  We  derive  the  needed  form  of 
the  Liouville  equation  following  a  procedure  suggested  by  Kozin  [5]. 

Consider  the  equations  of  motion  in  the  form  given  by  (1.27)  with  f 


x  *  Ax  , 


(1.27') 


where  x  is  the  (2nxl)  column  vector  whose  transpose  x1  has  the  form 

ip  <  • 

x  *  (q^ ,  ...,  q^;  q^,  ...,  q^}  and  A  is  the  (2nx2n)  matrix  given  by  the 

first  of  (1.28).  The  vector  x  is  the  state  space  form  for  representing 

the  system  response;  the  components  of  x  will  be  denoted  by  x^(t),  k  * 

1,  ...,  2n.  The  random  variables  in  A  are,  as  before,  denoted  by 

X. ,  ...,  X  .  We  will  rewrite  (1.27')  in  component  form  as 
x  m 


'  W#1 


mm 


m 
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*k  =  MV  X2n;  Xl*  Xm;  °  *  k  “  lf  **’  2n  *  (7,1) 

Let  p(x^ . x2n;  Xp  ...»  Xr;  t)  be  the  joint  probability 

distribution  of  the  random  quantities  x . .  x_  ;  X,,  ....  X  .  We 

1  2n  1  n 

define  the  characteristic  function  $  as 


*  =  E[exp  i(  E  9^x^(t)  +  E  ♦  X  )],  i 
k=l  j”l  3  J 


=  \|-1  . 


(7.2) 


The  differentiation  of  (112)  with  respect  to  time  gives 


U  -  E[i  E  0  x.  (t)exp  i(  E  0  x,(t)  +  E  ♦.X.)]  .  (7.3) 

k-1  *  k=l  K  K  j=l  J  J 

The  use  of  (6.1)  in  (6.3)  yields 

94  ^n  ^n  m 

-  i  E  9  E[g  exp  i(  E  0.x,  (t)  +  E  4>  .X .)]  .  (7. 

k=»l  K  X  k=l  K  K  j=l  J  J 

Since  (7.2)  is  essentially  the. Fourier  transform  of  the  joint  density 

function  p(x, ,  x„  ;  X,,  ...,  X  ;  t),  the  inverse  Fourier  transform  of 
l  zn  i  m 


(7.4) 


(7.4)  produces 


|E. 

j  8^ 


(7.5) 


The  solution  of  (7.5)  for  p  is  given  by  a  suitable  function  of  the 
independent  integrals  of  the  Lagrangian  system 


dt  -dp _  = 

1  =  3f,  3f?  fi 

p(3x+...+ax  ) 

1  2n 


(7.6) 


Let  u^,  ...,  u2n  be  2n-independent  integrals  of  (7.6).  Then  we  know 
that  the  general  solution  of  (7.5)  is 


P(V  X2n;  Xr  Xm;  °  “  h(V  U2u;  X1 . Xm;  °  L 

where  h  is  an  arbitrary  function  whose  form  is  determined  by  the  initial 


v  Vv.v'  .'V'/ 
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conditions  on  x. 

The  2n-integrals  u^ ,  U2n  may  be  determined  in  one  of  two  ways. 

First,  the  general  solution  of  (7.6)  is 


At 

x  =  e  x  : 

o 


-At 

e  x  =  u. 


(7.8) 


where  u  = 


(up  ...,  u2n) •  Further,  for  distinct  eigenvalues  of  A, 


-At  „-l 
e  =  U 


(7.9) 


where  X ^ ,  ...,  X^n  are  the  eigenvalues  of  A,  and  U  is  the  (2nx2n)  matrix 
of  corresponding  eigenvectors  of  A.  This  is  one  way  to  determine  the  2n 
independent  integrals  u  required  in  (7.7). 

The  second  method  employs  Laplace  transforms  and  the  equation 


x  *  Ax  with  x  =*  u  . 

o 

Then  from  the  Laplace  transform  of  this  equation,  we  find 


(7.10) 


I  ■> 

IT 

[  ■:> 


L  hds  -  A)  *]x  =  u  ; 


(7.11) 


it  should  be  noted  that  the  u's  depend  on  the  x's  and  the  X's.  To 
evaluate  the  inverse  Laplace  transform,  we  need  the  eigenvalues  of  -A 
or,  equivalently,  the  roots  of  (Is-A)  =  0. 


V.  s  -  sT  J  m*  f  f  \rj*  J  'j>  *  *  *  -  *»  1  »  *  r»  ”  -  '*  •  »  '  *  "  *  •  >  *  .>  • 
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Either  employing  (7.8)  or  (7.11),  we  have  the  required  integrals. 
Hence,  in  principal,  we  find  with  (7.7)  the  exact  expression  for  the 
joint  density  p,  the  arbitrary  function  h  being  determined  by  the 
specifics  of  the  problem  addressed. 

Consider  the  case  where  the  parameters  X^ ,  . . . ,  X  are  independent 
of  the  integrals  (i.e.  the  initial  vector).  Then  (6.7)  can  be  written 


P(X1,  x2n;  Xl*  Xm;  C)  “  hl(ul . u2n;  t)h2(Xl . V  * 

where  h^  is  the  joint  density  function  of  the  parameters.  Suppose  we 

are  interested  only  in  the  impulsive  response  for  q^;  this  means  at  t  =* 

0,  x  ■  1,  with  all  other  x.  *  0.  Then  at  t  *  0, 
n+1  j 


P(x1.  *  *  *  *  x2n’  Xl»  •••»  Xm;  0) 

*  <5(x1)...6(xn)<S(xn+1-l)...  -  -  5(x2n^h2^Xl*  ***’  Xm^  * 

In  view  of  the  fact  that  at  t  *  0,  u  *  x  ,  we  see  that 

o 

^(up  ...,  u2q;  t)  -  6(u1)...<$(un)<S(un+1-l)...6(u2n)  , 
where  6(.)  is  the  delta  function.  Finally, 


(7.12) 


p(x1,  *•*»  x2n;  Xl’  **•’  Xm;  fc)  (7.13) 

*  6(u1)...5(un)5(un+1-l)...<5(u2n)h2(X1,  ...»  X^)  . 

Let  us  now  illustrate  this  process  with  a  simple  example. 

Consider  again  the  undamped  one  degree  of  freedom  linear 

2 

oscillator.  Let  K/M  ■  ft  be  a  random  variable  with  sample  value  denoted 
2 

by  u>  .  The  first  order  system  representation  of  the  equation  of  motion 


'At 
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Equation  (7.5)  becomes 


X1  "  x2  "  fl  * 

*  2 

X2  "  x!  "  f2  * 


It  +  x2  ^  ‘  xin2  ^  "  0  • 


and  (7.6)  takes  the  form 


(7.14) 


(7.15) 


dt  dXl  _  _^2 

1  X2  x^2  ' 

We  now  find  the  integrals  of  (7.16);  they  are 


(7.16) 


Uj  *  Xjcos  Qt  -  jj-  sin  lit  , 

U2  =  xl^  s*n  ^  +  X2  COS  * 

2 

Assume,  for  example,  fl  has  a  discrete  distribution  given  by 


(7.17) 


h2(ui2)  -  E  Pi6(u>2  -  u>2)  . 


(7.18) 


If  we  are  interested  in  the  impulse  function,  at  t  *  0  we  have  x^  -  0, 


x2  =  1 .  There,  (7.13)  becomes 


p(x1,  x2,  id  ;  t) 
5(u1)6(u2-l)Zpi<5(u2-a)2) 


(7.19) 


2  2  2 
=  5(x,cos  ait - sin  (Dt)5(tDx.  sin  u>t+x~cos  u>t-l)£p  ,<5((d  -id.)  . 

1  (D  1  l  j  1  1 

Let  us  determine  the  mean  of  x^  to  illustrate  a  possible  use  for 


(7.19);  we  have 


2  2 

Efx^  -  /dx1/dx2/d(D  [x1p(x1,  x2>  <d  ;  t] 


Some  manipulation  (see  [5])  yields 


(7.20) 


;>■ 5>Siv  i:  S:  ::  v 
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E(,E{x1  jm,k,c}{Xj}  -  Q(t)  -  //x1P1dx1dx2 
EcE{x2|m,k,c}{x2}  -  m0^(t)  *  //x2Pidxidx2 
Differentiation  (partial)  of  these  equations  with  respect  to  time 


produces 


But,  from  (7.23'), 


mi,o "  ff*ijrdxidx2 

m0,l  "  ^X2  If  dxldx2 


8(8^)  3^2P2) 


(7.26) 


(7.27) 


The  substitution  of  (7.27)  into  (7.26),  the  employment  of  (7.22),  and 
integration  by  parts  of  the  resulting  terms  on  the  right  of  (7.26) 


finally  yields 


“1,0  "  m0,l  * 

k  c 

“0,1  "  "  m  ®10  "  m  m0,l  ’ 


(7.28) 


The  same  procedure  will  produce  the  equations  for  the  conditional 
2  2 

moments  E^Jx^ |m,k,c) ,  Ec{x2 |m,k,c} ,  etc.  We  note  that  for  the  first 
conditional  moments  we  could  have  taken  the  conditional  expectation  of 
(7.22)  to  produce  (7.28);  however,  this  procedure  only  applies  to  the 
first  moments. 

We  integrate  the  moment  equations  (7.28)  to  obtain  the  conditional 
moments  as  a  function  of  time.  On  multiplying  these  moments  by 
h2(m,k,c)  and  integrating  over  m,  k,  and  c,  we  obtain  the  moments  of  x^ 
and  x2« 

It  is  clear  from  the  above  discussion  that  the  Liouville  equation 
will  provide  the  exact  solution  for  the  joint  probability  density 
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function  p(x  ,  x0  ;  X.,  ...,  X  ;  t)  in  the  absence  of  external 

1  zu  1  m 

forces  provided  the  integrals  . u2n  can  be  obtained.  Further,  it 

provides  a  straight  forward  method  for  determining  the  moments  of  x  from 
which  means  and  variances  of  x  can  be  obtained. 

The  Liouville  equation  applies  when  there  are  no  external  forces. 

We  are  interested  in  the  case  when  external  forces  are  present,  of 
course.  Let  us  see  what  can  be  done  along  these  lines. 

The  Fokker-Planck  equation  is  the  natural  extension  of  the 
Liouville  equation.  This  equation  has  already  been  derived  above.  (See 
Section  2  and  also  [8],  [6b]).  We  confine  attention  to  the  case  in 
which  the  external  force  vector  f  can  be  obtained  by  passing  gaussian 
white  noise  through  a  stable  linear  damped  system.  We  have  as  equations 
of  motion,  conditional  on  M  -  m,  K  **  k,  and  C  -  c. 


dx1  =  x2dt  , 

dx.  *  -  —  xdt  -  —  x„dt  +  x_dt  , 

2  m  m2  3 

dx^  »  -gx^dt  +  dB  ,  Xj  3  0  at  t  *  0  . 


(7.29) 


where  we  have  employed  the  differential  notation  in  this  case,  set  f  * 
x^,  and  where  dB  is  the  Brownian  motion  increment  (see  Section  2)  with 


E{dB}  -  0  ,  E{ (dB)2}  -  a2dt  . 


(7.30) 


The  last  equation  of  (7.29)  represents  the  fact  that  the  excitation  is 

obtained  by  passing  a  gaussian  white  noise  through  a  linear  first  order 

stable  filter.  We  notice  that  for  the  Ito  system  (7.29) 

T 

x  *  {Xj ,  x2>  x^}  is  a  vector  Markoff  process  that  generates  a  Fokker- 
Planck  equation. 


imt 


It  can  be  shown  that  in  this  case  the  Fokker-Planck  equation  for 
the  conditional  probability  density  function  is 

-T7T-  m  ~  3^— (x.p.)  ~  -j— {(-  —  x,—  x_+^-  x-)p.}  ~  -J— (-6x,p,)  (7.31) 
3 1  dXj  2  1  8X2  m  1  m  2  m  3  rl  3x^  3  1 

+  £i!A 

2  Sx3 

We  observe  that  all  but  the  last  term  on  the  right  are  the  same  as  would 
have  occurred  in  the  Liouville  equation  in  the  absence  of  f.  Let  the 
conditional  moments  be 


kl  k2  k3. 


\.k2.k3  ‘  E(I1  x2  *3  1 
kl  kl  k3 

*  ///x1  *2  x3  PjCxj,  x2,  x3)dx1dx2dx3 
Then,  proceeding  as  in  the  development  of  (7.28),  we  find 


(7.32) 


ml,0,0  =  m0,l ,0 

m0,l,0  “  "  &  ml,0,0  ~  &  m0,l,0  +  &  m0,0,l  (7.33) 

m0,0,l  "  ~0mO,O,l 

On  multiplying  the  solutions  of  (7.33)  by  l^m.k.c)  and  integrating  out 
the  condition  in  these  three  conditional  moments,  we  finally  obtain  the 
moments  of  x  as  a  function  of  time.  We  obtain  in  analogous  fashion  the 
differential  equations  for  the  second  conditional  moments;  we  do  not  do 
this  as  the  steps  are  of  a  mechanical  nature  and  not  of  direct  interest. 
The  main  point  to  notice  is  that  differential  equations  for  the 


conditional  moments  of  x  can  be  obtained  when  an  external  force  is 


present  in  the  equations  of  motion  provided  this  force  is  produced  by 
passing  white  noise  through  a  suitable  filter. 
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It  is  important  to  point  out  that  for  any  gaussian  external 
excitation  the  solution  vector  is  gaussian  conditioned  on  the  random 
parameters.  Therefore,  all  conditional  moments  can  be  obtained  but  not 
as  easily  as  above  [79]  . 

The  Liouville  equation  enabled  us  to  obtain,  in  a  straight  forward 
manner,  the  exact  expression  for  the  conditional  probability  density 


function  p^.  Reference  to  (7.31)  suggests  that  it  will  be  much  more 


difficult  to  obtain  p^  from  this  equation  and  we  shall  not  pursue  this 


line  of  thought  further. 
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VIII.  Mean  Square  Approximate  Systems 

We  consider,  In  this  section,  a  technique  for  including  disorder  or 
parameter  uncertainty  that  follows  a  different  line  than  taken  in 
Sections  a)  and  b) .  Specifically,  mean  square  systems  are  employed.  We 
begin  by  introducing  this  type  of  system  [13,17], 

Let  us  begin  with  a  very  simple  example  in  which  there  is  no 


disorder  and  no  damping.  Let  the  coordinates  be  q  ,  ...,  q  .  Then, 

1  n 


with 


2T  ■  VjV  2v  ■  VjV  5U  ■  f3(t><,j  • 


(8.1) 


where  summation  convection.  Then  (1.24),  with  mass  coefficients 
included,  can  be  rewritten  as 


V  k  +  Vk  '  fJ(t) 


(8.2) 


Let,  with  f  constant, 


f^(t)  -  fojCos((Dt  +  $) 


(8.3) 


Then,  the  forced  motion 


q^  *  u^cos(o)t  +  <p) 


(8.4) 


satisfies 


(kjk  •  “  mjk)uk  ■  £oj  • 


(8.5) 


These  equations  state  that  given  the  f  and  to,  the  u^  are  determined  by 
the  solution  of  this  linear  system  of  equations.  Further,  if  to  is  the 


th 


natural  frequency  to^  and  the  u^  define  the  r  mode  shape  arjc»  then  the 


fQj  must  vanish.  Let  us  look  at  the  natural  frequency  problem  in  an 


unorthodox  manner. 


Suppose  we  pick  an  ui  and  a  set  of  u^  which  may  not  be  one  of  the 
natural  frequencies  and  normal  modes.  Then  the  right  of  (8.5)  will  not 
be  zero  and  we  need  force  amplitudes  to  produce  this  motion: 

(kjk  ‘  “Sk)uk  ‘  '  (8'6 

The  are  the  amplitudes  required  to  maintain  the  assumed  motion;  we 
regard  the  as  the  amplitudes  of  the  constraint  forces  required  to 
produce  the  motion. 

Consider  next 


I(n,u>)  «  Z  e*  >  0  .  (8.7 

1  J 

For  a  fixed  uj,  this  is  a  positive  definite  quadratic  function  of  the 

u's.  We  can  use  this  equation  to  find  the  natural  frequencies  and 

corresponding  normal  modes  ar^*  Assume  the  u's  are  normalized  in  some 

manner  (for  example,  u^  »  1  or  better  m^u^u^  -  1).  For  fixed  <d,  we 

find  the  minimum  of  I(u,  <u)  >  0.  Notice  that  if  u  ■  id  the  u's  that 

r 

2 

produce  a  minimum  are  the  a  and  I(a  ,  ,  w  )  =  0,  since  the  e  *  0,  1  « 

rk  rk  r  j 

1,  ...,  n,  in  this  case.  It  follows  that  if  for  a  specified  o)  we  find 

2 

the  minimum  of  I(u,  uj  )  and  this  minimum  equals  zero,  then  this  u  Is  a 

natural  frequency  and  the  u  that  produce  this  zero  minimum  are 

proportional  to  the  corresponding  normal  mode.  Let  us  consider  another 
interesting  aspect  of  this  method. 

Consider  a  frequency  window  g(u>)  with  the  following  properties: 


■k  u~m  „  >•  '  v  .  »  _  -  v,  -*'**►.•  ,  *  .  *  »  *  .  •  *-*-*-« 
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g(ui)  >  0  ,  u)  <  oj  <  a)  , 


(8.8) 


/  g(io)du)  -  1  ,  /  a)  g(aj)du  <  •  . 


Replace  (8.7)  with 


I(u,  g)  -  /  Ee.du)  . 

'  J 


(8.9) 


Find  the  u  that  makes  (8.9)  a  minimum.  The  interesting  feature  of  this 


method  is  that  if  there  is  a  natural  frequency  of  the  system  in  the 


frequency  interval  (<o  ,  m  ),  the  u  in  min  I(u,  g)  determine  the  normal 


mode  of  this  natural  frequency.  Let  these  u  be  in  component  form 


{u<r),  ...,  u^r^};  then  the  corresponding  natural  frequency  is 
in 


determined  by  the  Rayleigh  quotient: 


.(r)  (r) 


r  m  (r)  (r)  ’ 

mjkUj  "k 


(8.10) 


where  we  assume  we  have  found  the  r  normal  mode  and  its  natural 


frequency.  It  follows  that  if  there  is  concern  that  an  interval 


(<»>  ,  oj  )  contains  a  natural  frequency,  we  have  method  for  determining 


if  this  is  the  case  without  determining  all  natural  frequencies  of  the 


system.  References  [13,17]  give  details  on  this  matter  we  cannot 


discuss  in  this  Report. 


The  computational  problem  of  finding  the  minimum  of  I(u,  u>  )  is 


carried  out  using  one  of  a  number  of  computer  codes  based  upon  conjugate 


gradient  techniques,  and,  hence,  is  not  a  problem. 


So  far,  there  has  been  no  disorder  in  our  system;  i.e.  the 


parameters  m„  and  klt  have  been  assumed  to  take  definite  values.  Let 
jk  jk 
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us  assume  at  this  point  that  mass  and  stiffness  contain  random 


variables.  We  define,  in  this  case. 


I(u,g)  -  E{  /  g(aOZ^du>}  , 

'  1  2  ~ 


(8.11) 


where,  in  vector-matrix  form 


n  o  t  2  T  2 

=  u(K  -  w  M)  (K  -  m  M)u  . 

1  2 


(8.12) 


Since  E  only  operates  on  E £.  in  (7.11),  we  have 

1  2 


E{Ze2}  -  E{uT(K  -  m2M)T(K  -  w2M)u}  , 
1  J 


(8.13) 


and  a)  is  a  fixed  parameter  in  (8.13).  We  assume  the  u  are  parameters  to 


be  determined.  Thus,  (8.13)  takes  the  form 


E{Ee2}  -  uTE{(K  -  u2M)T(K  -  <j2M)}u  . 
1  2 


(8.14) 


2  T  2 

We  note  that  (K-o^M)  =K-iuM  because  of  the  symmetry  assumed  for  K 


and  M.  In  all  events,  means  and  second  moments  of  K  and  M  are  all  the 


information  needed  to  determine  the  expectation  in  (8.14). 


We  then  proceed  as  in  the  deterministic  case,  since  I(u,  g)  has  a 


deterministic  form. 


To  relate  (7.7)  to  (7.13),  all  we  have  to  do  is  assume 


g(m)  -  5(u  -  m)  , 


(8.15) 


where  5(.)  is  the  delta  function.  The  substitution  of  (8.15)  into 


(8.11)  yields 
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I(u,u>)  -  UTE{(K  -  u2M)T(K  -  0)2M)}u  .  (8.16) 

This  expression  differs  from  (8.7)  because  of  the  assumed  random 
parameters  in  K  and  M.  If  in  (8.7),  to  is  a  natural  frequency  of  the 
deterministic  system,  I(u,oj)  =  0.  The  I(u,w)  >  0  in  (8.16)  because  of 
the  random  parameters.  Use  of  this  fact  has  been  made  in  [32]  to  obtain 
an  estimate  of  the  variance  of  natural  frequency  u>^ ;  the  formula  is 

I(u,  w2) 

Var  u - ,  (8.17) 

r  4t/ 
r 

where  is  the  rC^  natural  frequency  and  uf  is  the  corresponding  normal 
mode  for  the  system  with  mean  parameter  values.  Monte  Carlo  simulation 
[32]  reveals  that  (8.17)  can  be  conservative  and  a  correction  is 
suggested.  Equation  (8.17)  is  easy  to  use  since  a  minimum  for  I  is  not 
required.  Further,  (8.17)  provides  a  much  simpler  method  for  estimating 
the  variance  of  than  given  in  Section  V.  However,  mean  square 
approximate  systems  do  not  provide  any  information  on  E{w^}  or  on 
variability  in  mode  shape.  Let  us  next  consider  how  these  systems  apply 
to  estimating  frequency  response  with  parameter  uncertainty  present. 

We  take  the  equations  of  motion  in  the  form 

•  •  • 

Mq  +  Cq  +  Kq  -  f  .  (8.17) 

The  frequency  response  Z  *(u>)  and  Z (w)  are  defined  by  (6.47)  and  (6.48). 
For  the  external  force, 

f  =  $jreiUt  *  (8.18) 

with  r  fixed  and  6  *  0  for  j  *  r,  <5^  =  1,  the  component  form  of 

(8.17)  is 


I 
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qj  =  Zjr(u,)e 


(8.19) 


which  is  exact. 


Suppose  we  try  to  approximate  (8.19)  with 


qj  ejr6 


(8.20) 


where  the  8  ,  are  not  known  in  advance.  The  equations  of  motion  now  are 
jr 


not  satisfied  and  we  must  introduce  constraint  forces  to  bring  about 


their  satisfaction  as  in  (6.1): 


(K.,  -  oi  M..  +  iu)C..)8,  -  6  -  e 

jk  jk  3k  kr  jr  jr 


(8.21) 


KB.  “>  =  E{E£jrejr}  , 


(8.22) 


where  asterisk  denotes  complex  conjugate.  This  I  is  just  like  (8.11) 


except  the  8  have  replaced  the  u.  We  find  the  8  that  make  (8.22)  a 


minimum,  denote  this  8  by  8.  Then,  8  *  {8^,  ...,  8^}  is  the  mean  square 


approximate  to  the  Z^(u)).  It  can  be  shown  that  if  the  system  is 


deterministic  (i.e.  contains  no  random  parameters)  the  8  are  exactly  the 


jr 


The  e  are  complex;  hence,  the  right  of  (8.22),  when  written  out, 


KV“Sl+1“CJl)Blr-V>  '  <8'23) 


It  follows  that  the  minimum  of  (8.22)  is  for  the  real  and  imaginary 


parts  of  8^r«  This  added  complication  poses  no  additional  computational 


problem  [39,  42] . 
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The  method  also  supplies  an  error  criterion  that  makes  it  possible 
to  judge  the  accuracy  of  the  6^. 

References  [39,42]  describe  in  some  detail  how  the  above  technique 
can  be  applied  to  estimating  the  frequency  response  in  a  number  of 
structures  with  specific  attention  being  paid  to  numerical  details  of 
the  computations.  Reference  [51]  also  describes  how  these  techniques 
can  be  applied  to  the  construction  of  a  sequence  of  approximants  for  a 
complex  system  by  starting  from  a  highly  constrained  initial  system  and 
gradually  relaxing  the  constraints.  In  these  three  reference,  extensive 
use  is  made  of  the  error  criterion  to  determine  when  the  estimated 
quantities  (usually  frequency  response)  are  sufficiently  accurate  for 
the  purpose  in  hand.  A  comment  on  what  mean  square  approximate  system 
provide  is  now  in  order. 

We  observe,  for  example,  that  these  systems  enable  us  to  estimate 
frequency  response  Zk*(u))  by  means  of  3kr(u>) .  The  8kr(w)  are 
deterministic  numbers  that  take  into  account  the  means  and  variances  of 
the  statistical  parameters  of  the  structure.  Thus,  the  8^(01)  provide  a 
deterministic  estimate  for  Zk^(u>).  In  the  form  given  above  and  in 
[39,41,42],  it  is  not  possible  to  obtain  statistical  information 
concerning  the  Z  ^(m).  However,  it  is  possible  to  employ  Monte  Carlo 
methods  to  obtain  estimates  for  the  8^(00)  given  the  parameters  are 
sample  values  to  obtain  sample  values  for  the  Bkr(a>)  from  which 
statistical  information  can  be  obtained. 

The  statistical  energy  approach  (SEA)  merits  mention  at  this  point 
since  it  also  employs  average  energy  concepts  [55,56,57,59,80], 
Basically,  SEA  estimates  the  average  flow  of  energy  from  one  part  of  a 


structure  to  another.  For  example,  if  there  is  energy  input  into  one 
part  of  a  complex  structure,  this  method  provides  an  estimate  of  how 
this  energy  flows  into  another  part  of  the  structure.  Thus,  it  is 
possible  to  estimate  average  vibrational  energy  present  in  any  part  of 
the  structure.  Information  of  this  type  is  frequently  the  only  type  of 
information  it  is  possible  to  obtain  about  the  response  in  an  extremely 
complex  structure  containing  a  large  number  of  undamped  natural 
frequencies  in  I  hertz.  In  so  far  as  we  know,  nothing  has  yet  been  done 
to  include  the  influence  of  statistical  parameters;  however,  the  work 
given  in  [59]  suggests  it  might  be  possible  to  do  this. 


IX.  Bound  Determination 

Sections  V  and  VI  discuss  techniques  for  response  estimation  in 
lumped  linear  systems  whose  parameters  take  time  Independent  uncertain 
values  in  a  probabilistic  setting;  thus,  it  is  possible  by  these 
techniques  to  obtain  statistical  information  about  a  natural  frequency 
and  its  normal  mode,  a  frequency  response,  and  a  response.  The 
techniques  discussed  in  Section  VII  provide  a  deterministic  frequency 
response  in  which  account  has  been  taken  of  the  statistical  properties 
of  the  parameters.  This  section  assumes  that  all  we  know  is  a  bound  on 
the  parameter  values;  we  then  will  be  interested  in  what  can  be  said 
about  bounds  on  the  quantities  of  interest. 

Needless  to  say  the  stability  and  control  of  systems  in  which  only 
bounds  are  known  on  parameter  -values  and  the  disturbances  continues-  to 
be  of  interest  to  those  working  in  automatic  control,  economic  analysis, 
and  stability  theory  [34,58,60,61,62].  Maximum  response  in  structural 
systems  has  also  been  and  continues  to  be  of  interest  to  structural 
engineers  [58,61,62];  however,  it  is  usually  assumed  the  parameters  in 
these  systems  are  known  exactly.  Let  us  consider  our  problem  from  the 
point  of  view  of  those  in  automatic  control. 

Let  us  consider  the  simplest  possible  linear  system  in  order  to  fix 
ideas.  Assume  the  system  is  described  by  the  first  order  linear 
differential  equation  in  the  single  variable  x: 

x  -  -(a  +  Aa)x  +  f(t)  ,  (9.1) 

where  f(t)  is  the  external  disturbance,  £  is  a  known  positive  constant, 
Aa  is  an  unknown  constant  with  -Aa  <  Aa  <  Aa,  Aa  is  a  known  positive 


constant,  and  a-Aa  >  0.  Let  the  initial  condition  on  x  be 


x(0)  ■  0  .  (9.2) 

We  are  interested  in  bounds  on  x  for  t  >  0. 

We  consider  as  nominal  solution  of  (9.1)  the  case  where  Aa  *  0; 
thus  let  x  satisfy 

• 

x  ■  -ax  +f(t)  (9.3) 

with  initial  condition  x(0)  “  0.  The  solution  of  (9.3)  is 

x  =*  /e‘a(t“T)f(T)dT  .  (9. A) 

o 

Introduce  the  error  function 

y  *  x  -  x  .  (9.5) 

Then  y  satisfies 

y  -  -(a  +  Aa)y  -  Aax(t)  ,  y(0)  =  0  .  (9.6) 

Define  a  Liapunov  function  V  by 

V(t)  -  |  y2(t)  .  (9.7) 

Then, 

V  “  yy  (9.8) 

*  [-(a  +  Aa)y  -  Aax]y 

2  — 

*  -(a  +  Aa)y  -  Aaxy  ; 

Let  Aa  *  -Aa.  We  shall  see  shortly  that  this  yields  an  upper  bound  on  V 
which  is  equivalent  to  bounding  the  errror  given  by  (9.5).  Then,  on 
taking  the  absolute  value  of  the  second  term  on  the  right  of  (9.8),  we 


obtain 


V  -  -ja-'Sajy  +  Aa|x(t)||y| 


which  from  (9.7)  gives 


V  <  -2(a  -  la)V  +  ‘Sa|x(t)  |  (2V) 


Next,  let 


Y\ 

it 


which  gives 


n  -  v 


We  now  can  write  (9.9)  as 


*  1  V 

‘  2  yl/2 


n  <  -(a  -  ~Sa)n  +  Aa  |x(t)| 


Let  us  solve  (9.11)  for  the  equality: 


n  «  -a n  +  -^  | x( t )  |  ,  n(0)  -  0 

\|2 


where  a  ■  a  -  Sa,  obtaining 


n  =*  ~  /e"-(t'T)|x(T)|dx  , 


\  2  o 


or,  by  (9.10)  and  (9.11), 


V1/2  <  /e"-(t_T)|x(T)|dT  (9. 

\  |  2  o 

Finally,  by  appealing  to  the  Schwarz  inequality  for  integrals,  we  can 


write  (9.14)  as 


_  .  1/2  1/2 
,1/2  <  fe“2a(t-T)dT 

\  1 2  o  o 


Let  us  examine  this  expression. 


]** 

si 


First,  the  first  integral  on  the  right  can  be  evaluated  as 


1  -2at. 

21  (1  -  e  ~  > 


1/2 


(9.16) 


Second,  the  left  equals 


iy(t) 

\|7 


•,  by  (9.7).  Hence  (9.15)  reduces  to 


|y(t) |  <  - — -  [1  -  e“2(a"Aa)t]1/2[/^2(T)dx]1/2  (9.17) 

\|2(a-Aa)  ° 

At  t  *  0,  we  get  |y(a)|  «  0,  as  we  should.  The  term  (9.16)  increases  to 
unity  as  time  increases.  The  term 


[/72(T)dt]1/2  (9.18) 

o 

_2 

will  increase  with  increasing  time;  if  x  (t)  approaches  0  as  t 
approaches  ",  then  (7.18)  will  increase  to  a  positive  limit.  The  bound 
is  proportional  to 

Aa 

_ ___  > 

\  1 2(a-Aa 

indicating  the  larger  Aa,  which  determines  the  bounds  on  a,  the  larger 
| y ( t ) | .  We  note  that  the  bound  on  |y(t)|  increases  with  time  using 
(9.16);  (9.14)  will  provide  a  smaller  bound,  of  course.  In  all  events, 
equation  (9.5)  makes  it  possible  to  assert  that  x(t)  lies  within  the 
bounds 


x(t)  ±  | y ( t ) |  .  (9.19) 

Thus,  the  technique  described  enables  us  to  put  bounds  on  the  response 
given  bounds  on  the  system  parameter. 

The  example  employed  above  is  for  a  first  order  linear  ordinary 
differential  equation.  References  {  ]  suggest  that  the  technique  can  be 


extended  to  vector  differential  equations  of  this  type.  However,  the 
details  will  have  to  be  worked  out  to  determine  if  this  promising 
technique  can  be  put  to  practical  use,  since  it  may  turn  out  that  the 
smallest  possible  bounds  are  too  large  to  be  of  practical  value. 
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X.  Physics 

The  physicists  have  for  a  long  time  been  concerned  with  the 
vibrational  properties  of  disordered  systems.  In  particular,  they  have 
been  interested  in  crystals  in  which  disorder  is  present  (See  [63-75]). 
One  type  of  disorder  is  called  "substitutional  disorder";  here  one  or 
more  atoms  in  a  regular  crystal  are  replaced  with  another  atom  or  atoms 
different  from  those  in  the  crystal  (while  the  organization  of  the 
crystal  is  not  changed).  The  other  type  of  disorder  is  called 
"topological  disorder",  in  this  case  the  basic  organization  of  the 
structure  is  changed.  A  moments  reflection  will  indicate  that  our  main 
concern  is  with  the  first  type  of  disorder,  where  member  properties  may 
change  but  not  the  structural  organization.  For  the  second  type  to 
occur,  structural  members  would  have  to  be  removed,  added  or  rearranged 
differently.  Much  progress  ha®  occurred  in  dealing  with  substitutional 
disorder;  much  more  modest  gains  have  been  made  for  topological 
disorder.  It  might  be  hoped  that  much  of  what  the  physicist  has  done 
could  be  adopted  In  toto  in  studying  the  response  in  substitutionally 
disordered  structures.  Based  on  a  relatively  short  study,  this  does  not 
appear  to  be  the  case  because  of  different  interests. 

We  are  interested  in  natural  frequencies,  normal  modes,  frequency 
response,  impulse  function,  etc.  in  systems  with  a  relatively  small 
number  of  degrees  of  freedom  as  a  rule.  By  and  large,  the  physicist  is 
interested  in  estimating  the  number  of  natural  frequencies  in  a 
specified  frequency  interval;  to  do  this,  a  frequency  spectrum  must  be 
estimated.  Hence,  they  are  interested  in  the  effect  of  various  degrees 
of  disorder  on  this  spectrum  from  which  optical,  thermodynamic. 
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electrical,  etc.  properties  are  obtained.  Thus,  our  Interests  are  very 
different  from  theirs. 

The  early  work  of  Hori  and  Asaki  which  introduced  the  method  of 
transfer  matrices  [68,69]  did  provide  a  technique  that  has  proved  useful 
when  dealing  with  chain  like  structures.  However,  mechanical  structures 
do  not  usually  possess  organizational  regularity  as  in  crystals.  Hence, 
as  noted  above,  transfer  matrices  have  a  limited  range  of  application. 

References  [73,74]  indicate  how  a  Green's  function  (impulse 
function)  can  be  employed  to  derive  information  about  the  frequency 
spectra.  So  far,  we  have  not  been  able  to  determine  how  this  technique 
might  apply  to  structures.  Nonetheless,  detailed  study  may  reveal  there 
are  possibilities  that  have  been  overlooked  in  this  brief  survey. 

The  early  work  of  Born  and  Brillouin  [3,6]  on  the  effect  of 
experimental  errors  on  the  motion  of  classical  mechanics  systems  is  of 
course  classic  and  is  worth  reading  just  for  cultural  reasons. 
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Appendix  A 


Zarghame's  method  with  mass  variability  Included. 

Assume  the  mass  matrix  M  can  be  written  in  the  form 
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where  M  is  the  mass  matrix  when  all  structural  numbers  take  on  their 
mean  masses,  and  M  the  mean  stiffness  matrix  of  the  element.  Now 
some  of  the  may  correspond  to  some  of  the  X^  in  (14).  Hence,  to 
include  this  possibility,  we  set 
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Equation  (19)  is  replaced  by 

•  • 
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where  q  is  the  same  (nxl)  column  vector  as  in  (20).  Equation  (23)  is 
replaced  by 

(K  -  u)2M)a  -  0  (5A) 

and  a  is  the  same  (nxl)  column  vector  as  in  (22). 

We  apply  the  same  procedure  as  employed  before  to  obtain  the 
partial  derivatives  of  and  with  respect  to  the  W^,  replacing  the 
orthogonality  conditions  (26)  by 
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